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1  Introduction 


Friction  and  rubbing  of  materials  are  among  the  most  common  phenomena  in  mechanics, 
occurring  whenever  two  solid  bodies  come  into  contact.  These  phenomena  are  responsible 
for  a  variety  of  occurrences  in  everyday  life.  Some  of  them,  such  as  tire  traction,  are  very 
useful;  others,  like  violin  music,  are  asthetic  and  pleasing,  and  many  others,  such  as  noises, 
vibrations,  and  wear,  are  extremely  unpleasant  and  deleterious  to  mechanical  systems.  This 
common  occurrence  of  friction  and  the  diversity  of  its  effects  underscore  the  extreme  impor¬ 
tance  of  a  deep  understanding,  and  the  need  for  modeling,  and  control  of  friction  phenomena. 
It  is  well  known,  however,  that  the  phenomena  of  contact  and  friction  of  solid  bodies  arc 
among  the  most  complex  and  difficult  to  model  of  all  mechanical  events,  primarily  due  to 
the  complex  structure  of  engineering  surfaces,  the  severe  elasto  plastic  delormat'.on.  damage, 
heat  generation,  atomic-range  interactions  that  take  place  on  typical  contact  surfaces,  the 
presence  of  contaminants,  lubrication,  and  even  chemical  reactions  on  t  hese  contact  surfaces. 

Efforts  toward  an  understanding  of  friction  phenomena  and  of  modeling  friction  began 
with  the  historical  works  of  Amontons  [2]  and  Coulomb  [22]  over  two  centuries  ago.  Since 
then,  an  extensive  body  of  experimental  and  theoretical  work  has  accumulated  on  general 
tribology,  and  a  good  empirical  understanding  of  the  subject  exists  today.  However,  t lie- 
progress  in  formulating  a  theoretical  background  and  designing  models  of  frictional  interfaces 
have  been  much  slower  to  evolve  than  experimental  investigations.  Although  considerable 
progress  in  this  direction  has  been  made  in  recent  years,  there  are  still  several  issues  that 
need  to  be  resolved  in  order  to  model  friction  and  predict  friction  phenomena  with  practical 
reliability.  One  of  the  most  difficult  problems  encountered  is  the  estimation  of  mate-rial 
constants  occurring  in  new  constitutive  models  of  frictional  interfaces.  These  difficulties 
reflect  an  urgent  need  for  constructing  new  constitutive  models  of  contact  and  friction  and 
for  estimating  the  necessarv  materia!  coefficients. 

Presently,  there  are  two  basic  approaches  for  the  development  of  mechanical  constitutive 
models  of  friction  and  two  resulting  types  of  frictional  interface  models.  These  are: 

1.  phenomenological  models  based  primarily  on  experimental  observations,  and 

2.  asperity- based  models,  formulated  via  a  theoretical  analysis  and  statistical  homoge¬ 
nization  of  the  microscale  deformation  of  surface  asperities  in  contact,  with  an  opposing 
surface. 


Unfortunately,  to  date,  none  of  these  approaches  has  produced  completely  satisfactory 
results.  It  is  well  known  that  experimental  results  depend  strongly  on  the  characteristics  of 
the  test  apparati,  so  the  results  of  different  tests  on  the  same  sample  can  be  considerably 
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scattered.  Moreover,  these  macroscopic  experimental  measurements  do  not  provide  sufficient 
insight  into  the  nature  of  the  phenomena  occurring  on  the  contacting  surfaces,  the  severity  of 
the  deformation,  propagation  of  damage,  etc.  Such  an  insight  can  be  provided  by  asperity- 
based  models,  which  are  based  on  the  microscale  analysis  of  deformation  and  the  relative 
sliding  of  surface  asperities.  However,  predictions  made  using  classical  or  existing  asperity- 
based  models  were  not.  generally  applicable  to  the  environments  normally  met  in  engineering 
applications.  The  main  reason  is  that  these  classical  models  were  based  on  analytical,  closed 
form  solutions  of  the  deformation  of  a  surface  asperity,  which  required  gross  simplifications 
of  the  geometry  of  the  asperity  and  of  the  constitutive  models  of  the  contacting  materials 
(elastic,  plastic,  or  at  most  elasto-plastic). 


1.1  Objectives  of  the  Project 

In  this  project,  a  new  approach  for  constructing  constitutive  models  of  friction  has  been 
developed  that  provides  a  realistic  link  between  microscale  phenomena  occurring  on  con¬ 
tacting  surfaces  and  macroscale  phenomenological  models  of  the  interface.  This  approach 
involves  the  use  of  special  finite  element  methods  in  the  modeling  of  complex  deformations 
of  asperities  of  arbitrary  shape,  with  realistic  nonlinear  constitutive  models  of  the  contacting 
materials.  The  technical  approach  for  the  evaluation  of  the  microasperity  based  models  of 
contact  and  friction  consists  of  two  stages: 

1.  Apply  the  finite  element  technique  to  analyze  the  nonlinear  mechanical  responses  of 
surface  asperities  of  different  heights,  shapes,  and  with  general  viscoelastoplast  ic  ma¬ 
terial  properties. 

2.  Apply  statistical  homogenization  techniques  to  evaluate  macroscopic,  phenomenologi¬ 
cal  constitutive  models  of  the  interface. 

The  approach  developed  here  provides  a  means  for  generating  a  variety  of  new  and  useful 
models  of  frictional  interfaces.  Depending  on  the  selected  level  of  complexity  of  the  model 
of  the  asperity,  a  viscoelastoplastic,  hyperelastic,  ur  brittle  material  can  be  considered,  the 
evolution  of  the  damage  to  the  surface  can  be  modeled,  and  the  effects  of  lubrication  and 
surface  contamination  can  be  taken  into  account. 

1.2  Research  Summary 

This  final  report  presents  the  results  of  research  work  performed  during  the  three- year 
project.  The  major  tasks  and  results  of  each  year  arc  briefly  summarized  in  this  section. 


The  first  year  of  the  project  was  dedicated  to  building  a  solid  foundation  for  the  statistical 
homogenization  procedures,  as  well  as  initial  work  on  the  finite  element,  modeling  capability. 
In  particular,  the  following  tasks  were  accomplished  in  the  first  year: 

1.  A  detailed  study  of  the  mechanics  and  statistics  of  asperity-based  models  of  contact, 
friction,  and  adhesion. 

2.  Formulation  of  a  complete  theoretical  background,  and  development  of  computer  codes 
for  the  calculation  of  macro-scale  interface  parameters  from  profilornetric  data  of  the 
surface  and  from  finite  element  analysis  of  a  family  of  representative  asperities. 

3.  Initial  work  on  the  research-type  finite  dement,  code  for  the  nonlinear  analysis  of  surface 
asperities  in  contact  with  the  opposing  surfaces. 

In  the  second  year  of  the  project,  the  work  focused  on  the  complete  development  of  the 
finite  element  asperity-modeling  code,  and  on  initial  tests  and  experimental  verification. 
This  included  the  following  tasks: 

4.  Development  of  the  three-dimensional  adaptive  finite  element  code  for  the  analysis  of 
surface  asperities  in  contact  with  opposing  surfaces.  The  starting  point  for  this  effort 
was  an  existing  in-house  finite  element  kernel,  which  was  extended  and  customized  to 
satisfy  the  objectives  of  this  project.  The  development  effort  focused  on  the  implemen¬ 
tation  of  elastic  and  viscoplastic  three-dimensional  solid  models,  on  the  development 
of  contact  and  sliding  resistance  algorithms,  as  well  as  on  extensions  of  the  graphics, 
user  interface,  and  adaptive  algorithms  needed  for  this  project. 

5.  Design,  development,  and  performance  of  Phase  I  of  the  verification  experiments,  which 
were  oriented  on  the  testing  of  numerical  models  of  nonelastic  surface  asperities  in 
contact  with  a  rigid  flat.  Special  custom-shaped  “asperities"  were  used  at  this  stage. 

6.  Verification  of  finite  element  asperity  models  by  comparison  with  the  analytical  and 
experimental  results.  The  numerical  predictions  were  compared  with  existing  analytical 
solutions  for  selected  simplified  cases  (Hertz  problem)  and  with  the  experimental  results 
obtained  for  fully  nonlinear,  elastoplastic  contact  problems. 

7.  Introductory  tests  of  the  complete  homogenization  procedure,  to  study  the  macroscopic 
behavior  of  homogenized  interfaces. 

The  third  and  final  year  of  the  project  was  dedicated  to  the  .actual  development,  of  new, 
asperity- based  constitutive  models  for  a  variety  of  interfaces,  and  to  their  comparisons  with 
analytical  and  experimental  results.  The  particular  tasks  completed  in  the  third  year  include: 
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8.  Evaluation  of  the  results  of  experimental  verifications  performed  in  the  second  year  of 
the  project. 

9.  Design  and  performance  of  the  Phase  II  verification  experiment.,  dedicated  to  detailed 
studies  of  the  behavior  of  real  engineering  surfaces  under  normal  and  tangential  loads. 

10.  Development  of  asperity- based  constitutive  models  for  the  surfaces  studied  in  the  ex¬ 
periment,  and  comparison  of  numerical  and  experimental  results. 

11.  Additional  studies  of  the  influence  of  surface  roughness  on  the  elasto-plastic  response 
of  the  interface,  as  well  as  introductory  modeling  ot  asperity  behavior  under  frictional 
loads. 

12.  Formulation  of  a  theoretical  background  for  the  application  of  asperity-based  models  of 
interfaces  in  the  modeling  of  dynamic  f fiction  phenomena.  This  included,  in  particular, 
development  of  analytical  formulas  (models)  representing  the  behavior  of  asperiiy- 
based  models. 

Additionally,  an  extensive  study  of  error  estimation  techniques  and  lip-adaptive  mesh 
refinement  strategies  was  performed  for  various  classes  of  problems  involving  contact  and 
friction. 

The  final  result  ot  this  project  is  a  proven  and  workable  approach  to  the  development  of 
asperity- based  constitutive  models  of  frictional  interfaces,  together  with  relevant  research- 
type  homogenization  software.  These  results  are  directly  applicable  in  the  analysis  of  a 
variety  of  friction  phenomena,  such  as  the  kinetic  coefficient  of  friction,  friction-induced 
noises  and  vibrations,  surface  compliance  for  bearing  applications,  real  contact  area  for 
electric  and  heat  interfaces,  and  introductory  studies  of  models  of  surface  damage  and  wear. 

1.3  Personnel 

I  he  research  effort  during  the  course  of  this  project  was  performed  by  a  highly  specialized 
team  of  COMCO  researchers.  The  principal  investigator  on  the  project  was  Dr.  .1.  Tinsley 
Oden.  President  and  Senior  Scientist  at  COMCO.  Assisting  extensively  on  the  project,  were 
Dr.  W.  VVoytek  Tworzydlo,  Director  of  Continuum  Mechanics  group,  and  Dr.  Witold  Cecot. 
Senior  Research  Engineer.  Additional  help  was  provided  by  Dr.  Jon  Bass,  Vice-President 
for  Research  and  Technology  and  Mr.  Olivier  Hardy.  Graduate  Research  Engineer. 

A  starting  point  for  the  finite  element  modeling  capability  was  a  proprietary  adaptive 
finite  element  kernel,  developed  by  COMCO  software  group. 

The  specialized  experimental  work  was  performed  by  Professor  C.  H.  Yew  of  the  Univer¬ 
sity  of  Texas  at.  Austin.  The  error  estimation  study  and  lip-adaptive  strategy  development 
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for  problems  with  contact  and  friction  was  performed  by  C.Y.  Lee,  Graduate  Student  at  the 
Texas  Institute  for  Computational  Mechanics  (TICOM),  University  of  Texas  at  Austin. 


1.4  Presentations,  Publications,  and  Technology  Transition 

1.4.1  Presentations 

The  research  related  to  new  models  of  contact  and  friction  was  presented  at  the  following 
professional  meetings: 

1.  113th  ASME  Winter  Annual  Meeting, 

Anaheim,  California,  November  S-13,  1992. 

2.  AFOSR  Grantees  and  Contractors  Meeting, 

“Research  in  Computational  Mechanics”, 

Washington  University,  St.  Louis,  May  20-21,  1993. 

3.  114th  ASME  Winter  Annual  Meeting. 

New  Orleans,  Louisiana,  November  28-December  3,  1993. 

1.4.2  Publications 

The  following  friction-related  papers  were  published  or  submitted  for  publication  during  the 
course  of  the  project: 

Ibrahim,  R.  A.  and  Soom,  A.,  Editors,  Friction-Induced  Vibration,  Cnatter,  Squeal,  and 
Chaos,  ASME,  De-Vol.  49,  New  York,  1992,  Tworzydlo,  W.  W.,  Becker,  E.  B.,  and  Oden. 
J.  T.,  “Numerical  Modeling  of  Friction-Induced  Vibrations  and  Dynamic  Instabilities",  pp. 
13-32. 

Wriggers,  P.  and  Wagner,  W.,  Editors,  Nonlinear  Computational  Mechanics  -  State  of  the 
Art,  Springer- Verlag,  Berlin,  1992,  Lee,  C.  Y.,  Oden,  J.  T.  ,  and  Ainsworth,  M.,  “Local  A 
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Tworzydlo,  W.  W.,  Becker,  E.  B.,  and  Oden,  J.  T.,  “Numerical  Modeling  of  Friction-Induced 
Vibrations  and  Dynamic  Instabilities”,  Applied  Mechanics  Reviews,  to  appear. 

Two  additional  papers  dedicated  to  asperity-based  models  of  contact  and  friction  a. re 
currently  in  preparation. 

1.4.3  Technology  Transition 

The  results  of  this  contract  and  previous  AFOSR-sponsored  contracts  dedicated  to  friction 
modeling  are  finding  their  way  into  practical  applications  in  engineering.  This  includes,  for 
example: 

Tire  modeling 

The  Oden-Martins  friction  model,  was  implemented  in  the  TIRE3D  tire  modeling  code, 
developed  by  COMCO  under  the  National  Tire  Modeling  Program  (NTMP).  The  code  is 
presently  being  used  by  NASA  and  Goodyear  for  analysis  and  design  of  rolling  tires. 

Modeling  and  Prediction  of  Friction-Induced  Noises 

The  results  of  AFOSR-sponsorecl  friction  projects  are  being  applied  in  practical  attempts 
to  understand,  model,  and  eliminate  friction-induced  noises  in  industrial  applications.  In 
particular,  Ford  Motor  Company  and  ORTECII  International  Research  Institute  are  using 
the  approach  developed  in  our  projects  to  eliminate  noises  in  automotive  components,  such 
as  the  squeaking  window  seal  in  the  Ford  Taurus.  Presently,  the  Computational  Mechanics 
Company  is  beuig  involved  in  this  team  to  provide  expertise  in  friction  modeling. 

Modeling  of  Earthquakes 

Recently  the  Computational  Mechanics  Company  was  awarded  a  research  grant  from  the 
U  S.  Geological  Survey,  for  a  project  dedicated  to  Modeling  and  Pn.diction  of  Earthquakes  as 
Unstable  Phenomena  of  Dynamic  Friction.  The  research  work  in  this  project  is  directly  based 
on  the  methodology  and  experience  developed  in  previous  and  present  contracts  sponsored 
by  the  AFOSR. 
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Modeling  of  Bearings 

Currently  several  potential  R  Sz  D  projects  are  being  discussed  with  major  bearing  man¬ 
ufacturers.  The  projects  under  consideration  will  apply  the  results  of  AFOSR-sponsored 
research  to  precise  modeling  of  contact  problems  inherent  in  bearing  design,  such  as: 

•  compliance  of  the  interface, 

•  real  contact  area, 

•  surface  wear  mechanisms,  etc. 

1.5  Outline  of  the  Report 

This  report  presents  the  results  of  the  second  year  of  effort  on  this  project ,  as  well  as  a  brief 
compilation  of  the  most  important  results  of  year  I.  In  particular,  Section  2  presents  a  study 
of  statistical  methods  of  homogenization  of  interface  parameters.  Of  particular  interest  are 
such  issues  as  extraction  of  surface  statistics  from  profilornetric  data,  calculation  of  an  asper¬ 
ity  distribution  for  random  surfaces,  and  the  practical  calculation  of  expected  macroscopic 
parameters  from  a  microasperitv  analysis.  In  Section  3,  a  detailed  formulation  of  the  bound¬ 
ary  value  problem  representing  the  deformation  of  a  surface  asperity  is  developed.  This 
formulation  includes  elastic  and  viscoelastoplastic  material  properties,  damage?  modeling,  a 
nonpenetration  condition  on  the  contact  plain?,  and  boundary  conditions  resulting  from  ad¬ 
hesion  forces  and  sliding  resistance  of  the  interface.  Section  4  presents  the  background  of 
the  adaptive  finite  element  technology  developed  for  the  analysis  of  the  deformation  of  a 
microasperity.  A  general  idea  of  the  hp- adaptive  finite  element  methodology  is  discussed  in 
this  section,  together  with  a  detailed  presentation  of  the  numerical  algorithms  used  for  the 
solution  of  elastic  and  viscoplastic  contact  problems. 

The  above  theoretical  part  of  the  report  is  followed  by  examples  and  tests  of  the  mi¬ 
croasperity  analysis.  In  particular.  Section  5  presents  some  basic  tests  of  numerical  models 
of  viscoplastic  material  behavior.  Then,  in  Section  6,  finite  element  models  of  asperity 
response  are  verified  by  comparison  with  the  Hertz  solution  and  with  experimental  mea¬ 
surements  performed  for  custom-shaped  asperities.  This  section  is  followed  by  studies  of 
asperity-based  interface  models  for  various  types  of  engineering  surfaces  (Section  7).  1  hen. 
in  Section  8,  detailed  comparisons  of  asperity- based  models  with  results  of  specially  de¬ 
signed  experiments  are  presented.  Following  Section  9  is  dedicated  to  studies  of  the  static 
coefficient  of  friction  for  the  interfaces.  Section  10  presents  studies  directed  towards  applica¬ 
tion  of  asperity-basd  interface  models  in  modeling  of  dynamic  friction  phenomena.  Finally, 
in  Section  11,  conclusions  of  this  work  are  summarized  together  with  remaining  research 
challenges. 


2  Asperity— Based  Models  of  Contact  and  Friction 


One  of  the  major  missions  in  tribology  is  the  development  of  constitutive  models  of  frictional 
interfaces.  Throughout  the  decades  a  variety  of  approaches  and  types  of  models  have  been 
developed.  They  can  be  classified  into  several  groups,  including: 

•  models  based  on  experimental  observations, 

•  microasperitv-based  models, 

•  phenomenological  models  developed  from  basic  principles  of  mechanics,  and 

•  models  of  the  type  related  to  plasticity  theory. 

It  should  be  noted  here  that  this  distinction  is  only  of  a  general  nature  and  most  of  the 
models  presented  in  the  literature  combine,  in  some  sense,  features  of  more  than  one  of  these 
groups. 

In  this  project  we  focus  on  the  development  of  new  asperity-  based  models  of  contact  and 
friction.  These  models  are  aimed  at  the  development  of  constitutive  equations  of  frictional 
interfaces  via  the  statistical  homogenization  of  the  deformation  of  surface  asperities  subject 
to  contact  with  an  opposing  surface.  The  advantage  of  the  asperity-based  models  is  that 
they  provide  good  quantitative  insight  into  the  phenomena  occurring  at  the  interface  and 
predict  additional  information  hardly  available  from  the  experiment,  based  laws,  such  as  the 
surface  plasticity  indices,  microfracture  indices,  etc. 

The  first  contact  model  that  was  constructed  to  predict  the  true  contact  area  can  be 
found  in  a  paper  by  Abbott  and  Firestone  [lj,  in  which  the  contact  surface  was  simulated  in 
a  network  of  spheres  that  are  truncated  upon  indentation  into  a  hard  Hat.  By  knowing  the 
hardness  of  the  softer  of  the  two  materials  in  contact,  an  estimate  of  the  true  contact  area 
could  be  made,  assuming  perfectly  plastic  deformations. 

An  important  advance  in  development  of  asperity- based  models  of  contact  is  represented 
by  the  pioneering  paper  of  Greenwood  and  Williamson  [44],  in  which  the  rough  surfaces  were 
viewed  as  a  randomly  distributed  population  of  elastic  asperities  with  randomly  distributed 
asperity  heights.  Each  asperity  was  assumed  to  be  spherical  and  elastic  and  its  deformation 
properties  governed  by  the  Hertz  solution  for  elastic  contact.  Experimental  evidence  was 
provided  to  support  the  assertion,  now  widely  held  by  tribologists,  that  for  normally  isotropic 
engineering  surfaces,  a  Gaussian  distribution  of  asperities  heights  generally  exists.  In  such 
models,  there  are  no  microfrictional  effects  on  the  asperities,  such  effects  leading  to  second- 
order  changes  in  contact  pressure,  a  result  established  nearly  two  decades  earlier  by  Mindlin 
[62].  In  a  related  paper,  Greenwood  and  Tripp  [13]  showed  that  contact  of  two  rough  surfaces 
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with  Gaussian  distributions  of  asperity  heights  on  which  asperity  contacts  were  misaligned 
was  equivalent  to  a  single  elastic  surface  with  a  Gaussian  distribution  of  asperity  heights 
impending  on  a  rigid  flat.  The  use  of  such  statistical  representations  of  surface  topography 
has  since  become  a  popular  approach  in  modeling  both  elastic  and  inelastic  contact. 

The  Greenwood-Williamson  model  was  based  on  the  assumption  that  only  the  asperity 
height  was  a  random  variable,  and  that  the  radius  R  of  each  peak  was  constant.  Several 
generalizations  of  such  random  topography  models  appeared  in  the  literature  of  the  1970s. 
The  paper  of  VVhitehouse  and  Archard  [95]  extends  the  random  -asperity  models  to  include 
random  heights  and  curvatures,  and  Nayak  [68]  provided  a  general  approach  to  random  sur¬ 
face  modeling  using  notions  of  joint  probability  distribution  functions.  In  t  his  same  vein,  we 
mention  the  work  of  Bush,  Gibson,  and  Thomas  [23],  who  derived  a  joint  probability  distri¬ 
bution  density  function  for  random  asperity  heights  and  curvatures  of  a  random  population 
of  elliptic  paraboloids  in  elastic  contact  with  a  smooth  rigid  flat. 

Such  random-microtopography  models  that  employ  a  deterministic  function  for  asperity 
peak  shapes  are  called  asperity  models.  One  source  of  possible  inconsistency  in  such  models 
has  to  do  with  the  fact  that  a  Gaussian  distribution  of  asperity  heights  and  curvatures 
for  a  given  asperity  shape  may  lead  to  a  non- Gaussian  cumulative  probability  distribution 
of  the  surface  height,  an  unrealistic  result  for  most  "engineering  surfaces.”  This  problem 
was  addressed  by  Hisakado  [49]  and  Hisakado  and  Tsukizoe  [50],  by  assuming  a  Gaussian 
PDF  (Probability  Density  Function)  for  surface  heights,  with  a  given  deterministic  asperity 
shape,  and  then  deriving  the  PDF  for  peak  heights.  Hisakado  [49]  assumed  a  paraboloidal 
asperity  shape  and  Hisakado  and  Tsukizoe  [50]  a  conical  shape.  Francis  [41]  points  out  that 
the  Hisakado  models  may  lead  to  unrealistic  PDFs  for  asperity  heights,  since  they  rnay  be 
strongly  dependent  on  the  asperity  shape  and  may  become  negative  for  paraboloidal  and 
conical  shapes. 

F.xtensions  of  asperity-based  models  to  microcontact  deformation  laws  involving  elasto- 
plastic  deformations  were  first  contributed  by  Hisakado  [49].  Hailing  and  Nuri  [47]  account 
for  plastic  deformation  of  the  interface  by  assuming  that  a  rough  surface  deforms  elastically 
while  contacting  a  nonlinearly  elastic  flat,  representing  strain-hardening,  with  each  micro¬ 
contact  defined  by  a  fully-plastic  spherical  indentation.  Significant  generalizations  of  these 
types  of  asperity  models  can  be  found  in  the  detailed  studies  of  Francis  [41],  who  intro¬ 
duces  the  notion  of  the  sum  surface ,  discussed  later  in  the  present  work.  This  enables  out' 
to  model  Gaussian  engineering  surfaces  with  asperity  shapes  that  a  paraboloidal  only  at 
their  vertices,  but  which  have  random  heights  and  curvatures,  using  the  joint  PDF'  of  Navak 
[68].  Moreover,  Francis  [41]  also  takes  into  account  elastic  and  fully  plastic  deformations, 
with  strain-hardening,  using  functions  determined  empirically  from  spherical  indentations 
of  various  metals.  We  also  mention  that  an  extension  of  the  Greenwood-Williamson  mode! 
of  spherical  asperities  with  Hertzian  elastic  contact,  constant  radii,  and  random  heights  to 
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cases  in  which  a  transition  to  perfectly  plastic  deformations  occur  was  recently  proposed  hv 
Chang,  Etsion,  and  Bogy  (26-28]. 

We  note  that  most  of  the  references  cited  above  dealt  with  attempts  to  model  either 
contact  without  sliding  motion,  or  purely  static  or  quasi-static  friction  effects. 

The  asperity-based  models  of  frictional  interfaces  are  constructed  in  five  basic  steps: 

1.  Perform  a  statistical  analysis  for  the  surface  profile  (profiles), 

2.  calculate  the  surface  statistics  (distribution  of  surface  height,  gradient,  and  curvature) 
from  one  or  more  set  of  profile  data, 

3.  calculate,  from  surface  statistics,  the  probability  distribution  and  density  of  surface 
asperities  of  different  heights  and  (possibly)  peak  shapes, 

4.  calculate,  by  analytical  or  numerical  methods,  responses  of  representatives  of  a  family 
of  surface  asperities  of  different  shapes  to  prescribed  load  programs, 

5.  calculate,  from  asperity  data  and  the  probabilistic  distribution  of  asperities,  the  ex¬ 
pected  values  of  the  interface  response  (normal  and  tangential  forces,  damage,  etc.)  to 
prescribed  load  programs.  This  response  characterizes  constitutive  properties  of  the 
interface. 

Several  variations  of  this  basic  scheme  may  be  derived  for  random  and  deterministic 
surfaces,  isotropic  or  anisotropic  finish,  etc.  In  this  case  a  general  classification  of  surfaces 
presented  (after  Nayak  [68])  in  Fig.  2.1  is  helpful. 

For  practical  purposes,  it  is  reasonable  to  consider  the  following  three  classes: 

(i)  Gaussian  isotropic  surfaces, 

(ii)  Gaussian  anisotropic  surfaces,  and 

(iii)  other  surfaces,  in  particular  deterministic  surfaces  obtained  by  special  finishing  tech¬ 
niques. 

The  flowchart  illustrating  the  homogenization  procedure  for  these  three  groups  is  pre¬ 
sented  in  Fig.  2.2. 

The  details  of  these  procedures  will  be  discussed  later.  Here  it  is  important  to  observe 
that  for  Gaussian  isotropic  surfaces  it  suffices  to  gather  profilometric  data  along  only  one 
profile  on  the  surface  and  to  consider  asperities  of  axisymmetric  peak  shapes.  For  Gaussian 
anisotropic  surfaces,  however,  one  needs  at  least  three  nonparallel  profiles  and  asperities  of 
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Figure  2.1:  A  general  classification  of  surfaces. 
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Figure  2.2:  Flowchart  for  statistical  homogenization  of  interfaces. 
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different  principal  curvatures  and  orientations.  Finally,  for  non-random  surfaces,  a  full  two- 
dimensional  map  2  =  z(x,y)  of  the  surface  may  be  needed,  and  asperities  may  have  various 
deterministic  shapes,  depending  on  the  surface  finish. 

2.1  Microstructure  of  the  Frictional  Interface 

We  begin  by  considering  the  contact  of  two  deformable  bodies,  I  and  II,  over  a  nominal 
contact  area  Ao,  as  illustrated  in  Fig.  2.3.  An  element  of  unit  nominal  contact  area  is  isolated 
for  study,  as  indicated  in  the  figure.  The  average  stress  vector  S  over  the  unit  contact  area 
has  components  of  force  P  and  Q  normal  and  tangential  to  the  unit  area,  respectively.  The 
situation  is  equivalent  to  that  of  two  typical  coupons  of  surface  material,  one  taken  from  the 
material  near  the  contact  surface  of  each  body,  pressed  together  with  a  force  P  normal  to 
the  tangent  plane  at  the  center  of  the  coupon  interface  and  simultaneously  subjected  to  a 
shear  force  Q  tangent  to  the  plane.  The  bulk  deformations  of  bodies  I  and  II  are  ignored, 
our  aim  being  only  to  characterize  the  mechanical  properties  of  the  contact  interface.  The 
nominal  unit  surfaces  in  contact  are,  for  the  present,  assumed  to  be  initially  flat  and  parallel 
to  one  another. 

It  is  standard  practice  to  depict  the  approximate  profile  of  rough  engineering  surfaces 
with  a  profilometer  or  stylus,  drawn  across  the  surface,  which  generally  yields  a  jagged  profile 
with  an  exaggerated  vertical  scale  of  the  type  shown  in  Fig.  2.4(a).  We  consider  two  such 
opposing  surfaces  1  and  2  which  are  to  ultimately  come  in  contact.  Refrence  planes  defining 
the  mean  asperity  height  of  each  surface  profile  are  established,  and  we  characterize  the 
shape  of  each  profile  by  introducing  functions  z\  and  z2,  given  the  height  of  asperities  above 
the  respective  reference  planes,  i.e.,  the  functions  2,-  =  2,(1,  y),  i  —  1,2,  with  (x,  y)  a  point  in 
the  parallel  mean-height  reference  planes,  define  the  profiles  of  the  rough  material  surfaces 
1  and  2,  respectively.  The  distance  h  between  planes  is  the  separation  of  the  surfaces,  and 
the  distance  between  actual  opposing  material  points  is  denoted  s.  Thus,  at  a  point  (x,y) 
on  the  reference  plane,  we  have 

s  ~  h  —  2  (2.1) 

where  2  is  the  sum  surface  (see  Francis  [41]), 

-  =  Ci  +  22 

Francis  Las  pointed  out  that,  from  the  fact  that  the  sum  2  of  the  surface  heights  appears 
in  the  geometric  relation  (2.1),  the  situation  is  equivalent  to  that  of  a  single  deformable 
surface  of  height  z  —  zi  +  z2  approaching  a  rigid  flat,  as  suggested  in  Fig.  2.4(b). 

Clearly,  the  undeformed  surfaces  overlap  whenever 

s{x.y)  <  0 
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As  the  normal  load  pressing  the  surfaces  together  increases,  the  separation  h  decreases 
and  at  each  minimum  of  the  function  s  a  microcontact  nucleates  and  expands  due  to  local 
deformation  of  the  surfaces. 

It  is  of  importance  to  note  that  arguments  presented  by  Francis  [41]  are  of  a  purely 
geometric  nature.  From  a  mechanical  point  of  view,  two  major  objections  can  be  raised 
here: 

1.  The  sum  of  two  asperities  (say,  spherical)  in  contact  with  a  rigid  flat  is  not  mechanically 
equivilant  to  two  spheres  in  contact— see  Fig.  2.5. 

In  particular,  the  distance-force  curves  P  =  P(a )  for  the  two  models  are  different. 
Moreover,  the  “asperity”  peak  on  the  sum  surface  corresponding  to  two  spheres  is  not 
spherical.  It  can  be  shown,  however,  using  the  Hertz  solution,  that  these  differences 
vanish  when  the  ratio  of  asperity  radius  R  to  the  contact  radius  r  goes  to  infinity 
(relatively  smooth  surfaces  at  moderate  loads). 

2.  In  the  case  of  contact  with  friction,  the  sum  surface  approach  will  not  model  the  friction 
component  due  to  the  interlocking  of  asperities.  Similarly  as  above,  the  importance  of 
this  effect  diminishes  with  increasing  surface  smoothness. 

In  view  of  these  remarks,  the  sum  surface  approach  seems  to  be  correct  and  justified  for 
typical  engineering  surface  finishes  at  moderate  loads.  Note  that  this  condition  is  also  re¬ 
quired  for  the  satisfaction  of  the  assumption  that  separate  microcontacts  do  not  interact 
mechanically  and  that  contacts  do  not  merge. 

It  is  well  known  in  tribology  that  techniques  used  to  produce  engineering  surfaces  usu¬ 
ally  produce  a  Gaussian  distribution  of  the  surface  heights  c,.  Moreover,  the  sum  z  of  two 
Gaussian  surfaces  is  also  Gaussian;  indeed,  Tallian  [85]  points  out  that  if  Z\  and  z2  are  not 
exactly  Gaussian,  their  sum  surface  will  be  closer  to  Gaussian  than  either  surface.  If  the 
shape  of  an  asperity  is  assumed  to  be  paraboloidal,  as  have  been  done  by  several  authors, 
then  the  peak  heights  and  curvatures  are  correlated  random  variables,  with  the  result  that  a 
Gaussian  distriution  of  heights  and  curvatures  may  lead  to  a  cumulative  probability  distri¬ 
bution  of  surface  heights  which  is  non-Gaussian.  This  issue  has  been  studied  by  Hisakado 
[49],  Hisakado  and  Tsukizoe  [50],  and  by  Francis  [41],  who  assert  that  if  the  peak  shape  is 
paraboloidal  only  at  its  vertex,  then  the  ensemble  of  peaks  can  be  made  to  conform  to  the 
Gaussian  distribution. 

2.2  Statistics  of  a  Random  Surface 

There  are  several  methods  of  homogenization  that  can  be  found  in  the  literature.  Few 
have  been  effectively  used  for  describing  nonlinear  frictional  phenomena.  As  one  possible 
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Figure  2.5:  Mechanical  difference  between  asperity  behavior  and  sum  surface  model. 
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technique  we  describe  a  general  approach  inspired  by  the  works  of  Lonquet-Higgins  [57-59], 
Nayak  [68],  and  Francis  [41];  see  also  Chang,  Etsion,  and  Bogy  [26-28].  To  fix  some  of 
these  ideas,  we  note  that  for  a  given  asperity  profile,  one  defines  the  autocorrelation  function 
C(X,  Y)  for  the  random  variable  z  =  z(x,y)  (the  surface  height), 

C(X,  Y)  =  lim  ~  f  f  z{x,y)z(x  +  X,y  +  Y)dx  dy 

a~*oo  ab  J o  Jo 
6— oo 

and  the  power  spectral  density  P(kx,ky)  as  its  Fourier  transform, 

]  roo  roo 

P(kx,  ky)  =  —  /  /  C(X,  Y)  exp  [-i(Xkx  +  Ykx))dX  dY 

47T"  J  —  oo  J  —  oo 

The  power  spectral  moments  are 

/-X)  roo 

/  P(kX,ky)k'xkldkX  dky  (2.2) 

-oo  J —  oo 

and  the  r.m.s.  roughness  a  is  the  variance, 

a2  =  m0 o  =  C(0,0)  =  f  [  P{kx .  ky)dkx  dky 

J  —oo  J  —  oo 

A  convenient  representation  of  a  continuous  random  surface  is  of  the  form  [58.68]: 

OO 

z{x,  y)  =  cos  (xkin  +  ykyn  +en)  (2.3) 

n=l 

where  amplitudes  Cn ,  wave  numbers  krn  and  kyn,  and  phase  zn  are  random  variables.  It 
is  assumed  that  there  are  an  infinite  number  of  wave  vectors  in  any  area  dkxdky  and  that  £ 
has  a  uniform  probability  density  in  the  range  (0,2tt).  The  power  spectral  density  is  related 
to  representation  (2.3)  by 

P(kX,ky)dkXdky  =  ~j2C2n 
“ 

the  summation  being  over  all  terms  with  (knx,kny)  lying  in  the  area  dkxdky  around 
( kx ,  ky).  The  power  spectral  moments  mpq  can  then  be  expressed  as 

=  (2--0 

~  n 
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Similar  definitions  and  representations  as  above  can  be  introduced  for  arbitrary  surface 
profile  2(3),  s  being  a  parameter  on  the  surface.  Of  particular  interest  are  spectral  moments 
of  a  profile  m0,  m2,  and  m4. 

2.3  Calculation  of  Surface  Statistics  From  Profile  Data 

Information  about  the  statistics  of  a  two-dimensional  surface  can  be  effectively  obtained 
from  profilometric  data  for  one  or  more  profiles  on  the  surface.  This  greatly  simplifies  the 
homogenization  procedure  because  both  experimental  measurements  and  statistical  post¬ 
processing  are  much  easier  for  one-dimensional  profiles.  In  this  section  we  discuss  details  of 
these  computations  for  both  isotropic  and  anisotropic  surfaces. 

2.3.1  Profiles  on  Gaussian  Isotropic  Surfaces 

It  was  shown  by  Lonquet- Higgins  [57-59]  and  Nayak  [68]  that  for  random  isotropic  surfaces 
the  mean  surface  height  and  non-zero  spectral  moments  are  expressed  in  terms  of  mean 
profile  height  and  profile  spectral  moments: 

~  surface  —  ~  profile 

m  00  -  m0 

m20  =  m02  -  m2 

3  m22  =  m04  =  m40  =  m4 

Therefore,  in  order  to  calculate  surface  statistics  it  suffices  to  perform  measurements  for  one 
profile  on  the  surface.  The  spectral  moments  of  the  profile  can  be  calculated  in  several  ways: 

•  from  the  definition  as  moments  of  power  spectral  density 

•  from  statistical  postprocessing  (sampling)  of  profile  data 

•  from  counting  zeros  and  extrema  of  the  profile 

Calculation  from  definition: 

m,  =  f  $(k)k'dk 

J  —  OO 

where  k  is  a  wave  number  and  ${k)  is  a  power  spectral  density.  This  particular  method  is 
rather  expensive  because  it  requires  evaluation  of  the  autocorrelation  function  and  the  power 
spectral  density  as  its  Fourier  transform. 
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It  is  much  easier  to  calculate  profile  spectral  moments  if  one  reinterprets  them  as  standard 
deviations  cr,a,<7  of  profile  heights  z,  slopes  z,  and  curvatures  z,  respectively: 

m0  =  cr2 
•  2 

m2  =  cr 

772.4  =  (72 


Assuming  that  the  profile  data  was  sampled  at  n  points  separated  by  the  interval  As  (see 
Fig.  2.6),  the  profile  statistics  can  be  calculated  from  the  following  sampling  formulas  [16]: 


(a)  mean  height,  slope  and  curvature: 


1  n 

5 

1  n 

5  ■  ;£* 
l  n 

?  ■  ;£*■ 


(b)  variations  of  height,  slope  and  curvature: 

1  n 
o  1  sr~ 


a2  - 


;2  _ 


n  —  1 

» —  i 

1  n 

=  *- 

U  1  t-1 

1  n 


n2 

J)2 

=)2 


The  values  of  first  and  second  derivatives  can  be  calculated  from  a  second  order  approx¬ 
imation  of  the  profile  shape,  discussed  in  Appendix  A.l.  Note  that  the  mean  slope  and 
mean  curvature  of  a  perfect  Gaussian  profile  should  be  zero.  For  real  profiles  they  may 
slightly  differ  from  zero.  Also  note  that  for  the  above  procedure,  shorter  wavelengths  can  be 
automatically  filtered  out  by  appropriate  selection  of  the  sampling  interval  As. 

An  ingenuous  alternative  way  of  calculating  m2  and  7n4  was  proposed  by  Lonquet-Higgins 
[58],  see  also  Nayak  [68].  The  densities  of  zeros  and  of  extrema  of  the  profile  are  expressed 
through  spectral  moments  as: 


20 


D 


zero 


-C^extr 


I  (—']  '2 

n  \m0J 

1  fm4\i 
7 r  \m2J 


1  & 

7T  a 

1  a 
tv  a 


By  zero  point  we  mean  here  the  point  where  2(3)  =  !.  The  counting  of  zeros  and  extrema 
can  be  performed  simultaneously  with  the  sampling  procedure  described  above.  Then,  after 
calculation  of  the  variation  <7,  variations  of  slopes  and  curvatures  can  be  obtained  as: 


CJ~  —  TV  (7  D  o 

a  =  TvaDcxir 


In  practice,  both  sampling  and  counting  methods  can  be  easily  implemented  in  the  same 
sampling  program.  Practical  comparisons  of  these  procedures  are  presented  further  in  this 
report. 

Another  parameter  necessary  for  the  homogenization  procedure  is  a  density  of  peaks  on 
the  surface,  defined  for  homogenous  surfaces  as: 


Dp  =  lim 

dz— co  d A 
dy—oo 

where  dA  —  dx  dy  is  the  surface  area  and  Np  is  the  number  of  asperity  peaks  within  this 
area.  The  density  of  surface  peaks  can  be  calculated  from  profile  parameters  [68]  as: 


D”  ~  6 *V3 


1  m4 


m  2 


2.3.2  Profiles  on  Gaussian  Anisotropic  Surfaces 

Basic  statistical  information  for  anisotropic  random  surfaces  consists  of  nine  moments  of  the 
power  spectral  density:  moo,  ■  •  • ,  m4o-  However,  since  the  properties  of  the  surface  do  not 
depend  on  the  orientation  of  the  x,y  axes,  only  certain  invariant  combinations  appear  in  the 
probability  distribution  of  the  surface  statistics  [59,68].  These  invariants  are: 


1.  moo 

2.  m02  +  m20 
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3.  m-2QmQ2  -  m2n 

4.  m40  -f  2m22  +  m04 

5.  m40mo4  -  4m13m31  -f  3m22 

6.  (m4Q  +  m22)  (m22  +  m04)  -  (m31  +  m^)2 

7.  -[m40  -  m?3)  -  m31  (m^m^  -  m13m22)  +  m22  (m31m13  -  m^2)  } 

From  three  profiles  in  three  nonparallel  directions  0;,i  =  1,2,3  nine  parameters  can  be 
defined:  m0(,),rn2(,),m4(,),z  =  1,2,3.  However,  since  m0(i)  =  m0(2)  =  m0(3),  then  these  three 
profiles  define  seven  constants — invariants  described  above.  This  means  that  three  nonpar¬ 
allel  profiles  suffice  to  define  surface  statistics  for  Gaussian  anisotropic  surfaces.  (Detailed 
equations  will  not  be  derived  here.) 

2.4  Calculation  of  Asperity  Statistics  From  Surface  Statistics 

The  primary  idea  of  asperity-based  interface  models  is  to  calculate  interface  parameters 
(normal  force,  friction  force,  etc.)  for  a  family  of  asperities  of  certain  deterministic  shapes  and 
to  obtain  expected  values  of  these  parameters  for  the  interface  from  a  statistical  distribution 
of  asperities.  This  requires  the  calculation  of  probability  density  of  surface  asperities.  For 
random  surfaces,  this  probability  density  can  be  expressed  in  terms  of  surface  statistics.  This 
problem  will  be  addressed  in  this  section. 


2.4.1  Asperity  Statistics  for  Gaussian  Isotropic  Surfaces 


For  Gaussian  isotropic  surfaces  two  random  variables  are  assumed  to  govern  the  distribution 
of  asperities:  asperity  peak  height  zp  and  mean  curvature  «.  The  joint  probability  density 
function  of  these  parameters  was  derived  by  Nayak  [68]  and  recast  in  a  different  form  by 
Francis  [41].  Here  we  present  the  formula  due  to  Francis: 


where 


/*«.■»  =  drr*  ~ 1  + 


T\/l  —  ft 


i  - 

£p 

a 

nondimensional  peak  height 

T]  = 

v/T5  * 

<7 

nondimensional  peak  curvature 

(3  = 

oa 

wavelength  spectrum  parameter 
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The  wavelength  spectrum  parameter  varies  for  random  surfaces  between  0  and  1:  zero 
corresponds  to  the  widest  wavelength  spectrum  (asperity  heights  and  wave  numbers  are 
not  correlated),  and  one  corresponds  to  the  narrow  spectrum  (longer  asperities  have  bigger 
heights).  Note  that  deterministic  surfaces  may  have  0  >  1,  see  Section  2.6.1. 

2.4.2  Asperity  Statistics  for  Gaussian  Anisotropic  Surfaces 

For  Gaussian  anisotropic  surfaces  the  representative  asperities  are  no  longer  of  axisymmetric 
shape.  Instead,  one  should  consider  asperities  with  elliptic  horizontal  cross  sections,  shown 
in  Fig.  2.7b. 

The  peaks  of  these  asperities  can  be  characterized  by  four  parameters: 

1.  zp  —  peak  height 

2.  «i,«2  —  principal  curvatures 

3.  a  —  orientation  of  the  main  axis  of  curvature 

Equivalently,  peak  height  zp  and  three  Cartesian  curvatures  (second  derivatives  of  z) 
kxx,  *xy,  Kyy  can  be  used.  The  joint  probability  density  function  based  on  all  these  param¬ 
eters  should  be  defined  as  f^mn2Q  (£,  */i,  <*)•  In  principle,  this  function  can  be  defined 

from  surface  statistics,  in  particular  spectral  moments  moo,  -  •  • ,  m^o  [59,68].  To  the  author's 
knowledge,  no  such  formula  is  presently  available  in  a  closed  form. 

2.4.3  Asperity  Statistics  for  Deterministic  Surfaces 

Some  special  types  of  finish  may  produce  surfaces  of  non- Gaussian  random  distribution  or 
deterministic  distribution.  For  such  arbitrary  surfaces  the  distribution  of  asperity  peaks 
cannot  be  obtained  from  profile  data  and  need  to  be  calculated  directly  from  a  surface  map. 
In  this  section  we  present  a  simple  sampling  procedure  to  calculate  asperity  statistics  from 
surface  data. 

We  assume  that: 

•  The  surface  is  homogenous. 

•  The  function  z(x,y)  (surface  height)  for  the  surface  is  given.  This  can  be  obtained 
from  two-dimensional  sampling,  holography  or  other  methods. 

•  Asperity  peaks  are  characterized  by  the  peak  height  zp  and  mean  curvature  k. 
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a)  b) 


Figure  2.7:  Typical  asperity:  (a)  isotropic  random  surface,  (b)  anisotropic  random  surface. 


25 


The  probability  density  of  asperity  height  and  curvature  f2pK(zp ,  k)  can  be  obtained  from  the 
following  sampling  procedure:  Cover  the  domain  SI  with  a  regular  mesh  of  points  (x,,  y,),  i  = 
l,n  presented  in  Fig.  2.8.  The  mesh  spacing  h  can  be  defined  so  as  to  filter  out  high- 
frequency  noise.  The  sampling  is  performed  by  looping  through  the  points  (i,,  y,),  i  =  l,n 
and  for  each  point: 

1.  Check  if  the  point  is  a  peak  or  near  a  peak  within  resolution  h.  The  peak  is  identified 
as  a  peak  if  its  height  z(xp,  yp )  is  greater  than  all  its  nearest  neighbors  (eight  for  interior 
points).  Alternatively,  more  elaborate  criteria  may  be  used. 

2.  If  the  point  is  a  peak,  then  calculate  the  second  derivatives  zIX,zyy  and  :xy.  Here 
simple  finite  difference  formulas  may  be  used  or  a  generalized  minimization  procedure 
as  presented  in  Appendix  A. 


The  mean  surface  height  and  standard  deviation  of  the  surface  height  are  calculated  as: 


The  joint  probability  density  of  asperity  peak  heights  and  curvatures  can  be  calcu¬ 
lated  after  locating  all  the  peaks  by  dividing  the  range  of  peak  heights  and  curvatures 
[•Zpmim  ’pm»*]  *  [Kmini  Kmax]  into  area  elements  AzpAk  (see  Fig.  2.9).  Then  for  each  area 
element  with  a  center  point  (cpi,  kj)  define 

f:K{~pi1Kj)  =  j) 

Jyp 

Here  Np  is  the  total  number  of  peaks  and  np(i,j)  is  the  number  of  peaks  within  the  area 
element  AzpAk,  identified  by: 


~pmin  “b  ( 1  1)A~  ^  ~p  -pniin  d”  lA~ 

^•min  (j  1  )Ak  ^  Kp  <C  S.'niin  "b  J  Ak 


The  values  above  define  the  discrete  values  of  the  joint  probability  density  function  of  as¬ 
perity  peak  heights  and  curvatures.  This  function  may  then  be  regularized  by  an  application 
of  appropriate  approximation  techniques.  A  similar  procedure  can  be  used  if  one  chooses 
to  characterize  asperity  peaks  with  more  than  two  parmeters.  such  as  cp,K],k2  and  a  for 
anisotropic  surfaces. 
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Figure  2.8:  Sampling  of  arbitrary  surface. 


K,j 


Figure  2.9:  Accumulation  of  joint  probability  density  fSK(zp,K)  for  arbitrary  surfaces. 
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2.5  Calculation  of  Macrocontact  Expectations  of  Interface  Pa¬ 
rameters 

Once  we  know  the  probability  distribution  of  asperity  heights  and  shapes  as  well  as  the 
values  of  any  parameter  X  for  single  asperity,  it  is  possible  to  calculate  the  expected  value 
of  X  for  the  interface.  This  procedure  varies  somewhat  depending  on  the  classification  of 
surface  type.  In  this  section  we  present  a  detailed  procedure  of  the  expectation  calculation 
for  Gaussian  isotropic  surfaces  and  outline  its  extensions  to  other  surface  types. 

2.5.1  Expectation  Calculation  for  Gaussian  Isotropic  Surfaces 

For  Gaussian  isotropic  surfaces,  the  following  parameters  are  needed  to  calculate  the  expected 
value  of  the  macroscopic  interface  parameter  X : 

i)  Np  The  number  of  peaks  within  the  contact  area  /10- 

ii)  fzv{£,  ??)  Probability  density  of  asperity  peaks  of  (nondimensional)  height 

£  and  mean  curvature  r?.  By  a  simple  change  of  variables  one 
can  define  f^v(zp,K). 

iii)  X(zp,  K,a,d)  The  value  of  parameter  X  for  different  peak  heights  zp  and  cur¬ 

vatures  k,  subjected  to  a  normal  approach  a  and  sliding  distance 
d. 

iv)  cr,  d,cr  Deviations  of  profile  heights,  slopes,  and  curvatures  used  to 

nondimensionalize  peak  heights  and  curvatures. 

The  expected  value  of  X  per  asperity  is  calculated  as 


E[X(a,d)'j  =  J_^J0  X(zr,K,a,s)f(„(zp,K)dKdzf 


(2.5) 


and  the  macrocontact  expectation  of  X  is 

X(a,  s)  =  NpE(X(a,d)) 

Note  that  even  for  Gaussian  isotropic  surfaces  there  exist  asperities  of  non-axisymmetric 
cross  sections.  However,  due  to  isotropy,  it  suffices  to  consider  only  axisymmetric  represen¬ 
tatives  of  certain  mean  peak  curvature  k.  In  this  project  we  choose  the  asperity  to  be  a 
cosine  hill  defined  in  a  local  coordinate  system  as: 

z(i,  y)  =  C  cos  kx  cos  ky 
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For  this  asperity  the  peak  height  and  mean  curvature  are  defined  as 

zp  =  C 
k  =  Ck2 

The  above  shape  is  consistent  with  a  generic  representation  of  a  Gaussian  surface  presented  in 
formula  (2.3).  This  is  different  than  approaches  presented  to  date  in  the  literature,  in  which 
asperity  peaks  were  usually  assumed  to  be  spherical  or  paraboloidal.  This  was  because  these 
works  were  based  on  analytical  solutions  for  asperity  deformation,  such  as  Hertz’  solution. 
In  this  work  we  are  modeling  the  asperity  by  the  finite  element  method,  so  it  is  possible  to 
use  the  model  which  does  not  suffer  from  inconsistencies  of  spherical  or  paraboloidal  asperity 
peaks. 

Although  for  Gaussian  isotropic  surfaces  the  probability  density  of  asperity  heights  and 
curvatures  is  analytic,  the  values  of  X(zp,rc, . . .)  that  we  obtain  from  finite  element  com¬ 
putations  are  not.  Therefore  the  expected  value  of  X  must  be  calculated  using  numerical 
quadrature.  This  quadrature  was  implemented  under  the  following  assumptions: 

(a)  The  domain  of  integration  is  truncated  to  the  subregion  [~min,  ‘max]  x  [«min,  *imax]< 
where  the  probability  density  is  large  enough  to  effectively  contribute  to  the 
final  integral  E(X).  This  region  is  defined  adaptively  (see  Section  2.6.2). 

(b)  The  parameter  values  X(zp, «,...)  are  given  (calculated  by  FE\1)  for  certain 
selected  values  of  peak  heights  and  curvatures  in  the  domain  of  integration.  These 
points  are  not  necessarily  regularly  distributed  within  the  domain  of  integration. 

The  numerical  procedure  for  the  calculation  of  the  integral  consists  of  the  following  steps: 

1.  Divide  area  {zmin,Zmax]  x  {Kmin,«max]  into  area  elements  A z  x  A k  (see  Fig.  2.10). 

2.  Calculate  the  integral  by  looping  over  cells  and  applying  numerical  quadrature  (trape¬ 
zoidal,  Simpson,  Gauss,  or  any  other)  according  to  the  formula 


JV 


E(X)  =  £ 


1=1 


m 

y  i  ■  ■  ■)  fti  ("p(°)’  K(a ))  a*  A«| 


Lor  =  l 


(2.6) 


where  i  is  the  number  of  integration  cells,  a  is  the  number  of  quadrature  points  within 
a  cell,  and  wa  is  the  corresponding  weight  factor. 
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Figure  2.10:  Numerical  integration  of  expectation  value  E(X). 
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Note  that  this  integration  requires  the  value  of  X  at  integration  points  (zp(a),  K(c))  within 
each  cell.  Since  it  may  be  difficult  to  perform  finite  element  analysis  for  values  of  zp  and 
k  corresponding  exactly  to  all  the  integration  points,  these  values  are  calculated  from  the 
original  data  points  using  the  error  minimization  procedure  presented  in  Appendix  A.  The 
quadrature  rule  currently  implemented  for  integration  are  the  trapezoidal  and  four-point 
Gauss  rule. 

Note  that  the  above  procedure  introduces  error  due  to  truncation  of  the  integration 
domain  and  due  to  numerical  integration.  This  leads  to  a  rather  unwelcome  result  that  even 
for  constant  X ,  the  calculated  expected  value  E(X)  would  be  different  than  X .  In  order  to 
compensate  for  this  error,  we  additionally  calculate  the  integral  of  the  probability  density 
(which  should  be  one): 


jV 


t=l 


m 

T.  fin  ( *■(<»)»  ^p(or))  WqXzXk 
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Then  the  corrected  value  of  expectation  of  X  is  calculated  as: 


E(X)  =  E(X)/I 


(2.7) 


This  procedure  assures  that  for  constant  X  the  expected  value  E( X)  is  equal  to  X. 


2.5.2  Expectation  Calculation  for  Random  Anisotropic  Surfaces 


As  mentioned  previously,  for  anisotropic  Gaussian  surfaces  one  has  to  consider  asperities  of 
random  peak  heights  zp,  principal  curvatures  «j  and  «2,  and  orientations  of  the  principal 
axis  a.  The  calculation  of  expected  values  of  the  interface  parameters  is  similar  to  equation 
(2.5): 


E 


K  i ,  K2 }  CL 
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(2.8) 


Note  that  presently  there  exist  no  closed  form  solution  for  the  probability  density 
of  asperity  peaks  of  random  heights,  principal  curvatures  and  orientations.  Also  note  that 
in  order  to  span  the  integration  space,  one  would  need  to  obtain  finite  element  solutions  for 
a  large  family  of  asperities,  corresponding  to  various  combinations  of  zp,  Ki,  and  o.  This 
would  be  a  very  expensive  task  computationally,  and  will  not  be  considered  in  this  project. 
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2.5.3  Expectation  Calculation  for  Deterministic  Surfaces 

The  calculation  of  expectation  values  E(X )  for  non-random  surfaces  follows  essentially  the 
same  numerical  procedure  as  for  Gaussian  isotropic  or  anisotropic  surfaces.  Depending  on 
the  surface  type,  the  peak  height  zp,  curvature  k  and  other  parameters  may  be  selected  to 
represent  typical  asperities.  The  joint  probability  density  „.(zp,  k,  . . .)  can  be  obtained 
from  surface  sampling  as  discussed  in  Section  2.4.3. 

2.6  Numerical  Verification  of  Statistical  Postprocessing 

The  homogenization  procedures  discussed  in  the  previous  section  for  Gaussian  isotropic 
surfaces  were  used  as  the  basis  for  the  implementation  of  specialized  software  for  this  purpose. 
In  this  section,  certain  basic  tests  of  this  software  will  be  presented. 

2.6.1  Verification  of  Profile  Postprocessing 

The  program  for  profile  postprocessing  was  designed  to  read  in  the  data  c(s)  for  one  or  more 
profiles  on  the  surface,  and  use  them  to  calculate  mean  profile  height,  slope  and  curvature 
as  well  as  spectral  moments  and  peak  density  for  the  surface.  Both  statistical  sampling  and 
counting  methods  were  implemented  (see  Section  2.3.1). 

Example  1 

In  the  first  example  a  deterministic  cosine  profile  was  generated: 


r(s)  =  C  cos  ks 

with  C  =  1  and  k  =  2.  For  such  a  profile,  the  mean  profile  height,  slope  and  curvature 
are  zero.  The  standard  deviations  of  height,  slope  and  surface  curvature  can  be  calculated 
exactly  to  be: 


a 


a 


a 


-7=  =  0.7071 

C 

k—=  =  1.4142 

k2-^=  =  2.8284 

v/2 
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The  corresponding  wavelength  spectrum  parameter  is  3  =  1.2247. 

The  above  profile  was  sampled  with  the  interval  0.2  which  corresponds  to  about  100 
points  per  one  wavelength.  The  results  of  this  sampling  are  presented  in  Fig.  2.11. 

The  calculation  of  mean  surface  heights,  slopes  and  curvatures  (which  should  all  be  zero) 
is  very  accurate.  So  is  the  calculation  of  the  standard  deviation  of  surface  height.  Deviations 
of  surface  slopes  and  curvatures  are  less  accurate  (up  to  7  percent  error),  which  is  caused  by 
an  approximate  calculation  of  slopes  and  curvatures.  For  this  particular  case,  the  counting 
method  gives  better  results  than  the  statistical  sampling  method. 

Example  2 

In  the  second  example  we  have  generated  a  quasi-random  profile  by  using  a  one-dimensional 
version  of  the  formula  (2.3).  A  series  with  40  components  of  a  quasi-random  distribution  of 
C„k,  and  £,  were  specified.  The  resulting  profile  is  presented  graphically  in  Fig.  2.12a. 

For  a  fully  infinite  series  and  random  C\,k,  and  c,  the  power  spectral  moments  are 
expressed  by  a  one-dimensional  equivalent  of  formula  (2.4).  For  truncated  series  this  is  not 
true,  but  a  reasonable  approximation  can  be  expected  for  the  values  of  spectral  moments. 
These  values,  calculated  for  the  above  profile,  are  shown  in  Fig.  2.12b.  The  profile  was  then 
sampled  using  the  procedure  described  in  Section  2.3.1.  The  results  are  shown  in  Fig.  2.12c. 

ft  is  somewhat  more  difficult  to  verify  the  results  in  this  case,  since  neither  of  the  methods 
produce  exact  results.  It  can  be  noted,  though,  that  the  results  for  surface  height  (mean, 
height  and  deviation)  are  the  most  accurate,  while  for  curvatures  the  differences  reach  up 
to  about  20  percent.  Other  tests,  not  discussed  here,  show  that  with  a  decreasing  sampling 
interval,  these  differences  become  smaller  (but  they  do  not  vanish). 

2.6.2  Verification  of  Expectation  Calculation 

To  perform  a  basic  test  for  the  calculation  of  joint  probability  density  of  surface  asperities 
and  for  expectation  values  of  A',  we  considered  the  following  set  of  test  data: 

1.  The  profile  and  surface  statistics  were  taken  from  Example  2  in  the  previous  section, 
in  particular: 
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sampling  results  for  single  wave,  C=l,k=2. 

+++T+++++++  ++++++++4~»--r +++++  T+rr++T++r+  +  +  +  J- 


Sampling  interval:  .2 

RESULTS 


Method  1 
( sampling) 


Method  2 
(counting) 


Mean  height :  -2 . 
Mean  siope:  -3. 
Mean  curvature:  -4. 

St.D.  of  height:  7. 
St.D.  of  slope:  1. 
St.D.  of  curvature:  2. 
Surface  peak  density:  1. 
Wavelength  spectrum:  1. 


3846e— 35 
343ie-05 
3874e-04 

07l5e-Jl 

3629e+00 

1 . 4143et00 

S272e+00 

2 . 7222e+00 

138  4e-0 1 

1 .2222e-01 

2243er00 

1.2727e+00 

Figure  2.11:  Sampling  result  for  a  simple  cosinusoidal  profile. 
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Sampling  interval:  .005 


'ACS  STATISTICS  r?.CM  SERIES 


Wavelength  spectrum.: 


1 . ISSOe-OO 
5  . 5221S-01 
4 . 33  ire-01 
7 . 4954e-01 


b) 


E  5  j 


Mechod  1 

Me choc  2 

(sampling) 

(countrr.g) 

Mean  height : 

S.7472e-02 

Mean  slope: 

-2 . 49 lle-03 

Mean  curvature : 

-i.3438e-02 

C) 

St.D.  of  height: 

1.12S9e+00 

St.S.  of  slope: 

6 . 5342e-0 1 

7 . 0748e-01 

St. D.  of  curvature: 

5  .  9455e-01 

5 . 1324e-01 

Surface  peak  density: 

2 . 535  9e-02 

I . 3898e-02 

Wavelength  spectrum: 

7. 3117e-01 

1.0  60  9e-*-00 

Figure  2.12:  Profile  sampling  test:  (a)  quasi-random  surface  profile,  (b)  surface  statistics 
from  series  expansion,  and  (c)  surface  statistics  from  sampling. 
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a  =  1.1620 

a  =  0.59221 

a  =  0.49315 

Np  =  0.025359 
Ao  =  1 


2.  The  area  of  integration  of  expectation  values  E(X)  was  defined  by: 


‘■pinin  5  ,  -pmax 

^min  —  0  ,  Kraax 


10 

5 


This  domain  was  subdivided  into  20  x  20  integration  cells,  with  trapezoidal  integration 
within  each  cell. 

3.  The  variable  X  was  assumed  to  be  identically  equal  to  one  (so  that  the  expectation 
value  should  be  one).  This  was  implemented  by  specifying  eight  data  points  with  a 
value  1.0,  randomly  distributed  within  the  integration  domain. 


The  results  of  the  numerical  calculation  of  expectation  values  are  presented  in  Fig.  2.13. 

Note  that  the  program  estimates  the  effective  support  of  the  probability  density  function 
ftv(zp,n),  which  is  identified  as  the  loci  of  points  where  fir)(zpK)  is  greater  than  10-4.  This 
is  done  to  avoid  integration  over  too  large  a  domain.  The  above  estimate  is  still  very 
safe — for  example,  for  the  profile  considered  here  the  estimated  effective  support  of  f(v{zpK) 
corresponds  to  peak  heights  between: 


^p  min  —  3-3  Zp  max  —  7.0 

while  the  real  profile  had  a  maximum  peak  height  of  only  2.5  and  a  minimum  peak  height 
of  about  —0.8. 

The  value  of  expectation  E(X)  calculated  by  integration  (formula  (2.6))  is  1.0043.  This 
is  an  effect  of  truncation  of  the  integration  domain  and  of  numerical  integration.  After 
correction  according  to  formula  (2.7),  the  value  of  expectation  E(X)  is  exactly  one. 

This  simple  example  confirms  the  correctness  of  the  theoretical  formulation  and  the 
software  used  for  calculation  of  joint  probability  density  ftv{£,ri)  and  for  the  expectation 
value  of  X  for  Gaussian  isotropic  surfaces. 
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INPUT  ECHO 

zpmin=-5 . OOOOe+OO  zpmax=  1.0000e+0l  kapmin=  0.0000e+00  kapmax=  l.OOOOe-Ol 
nksi=20  neta=20  irule=l 

sigma=  1.1620e+00  sigdot=  5.9221e-01  sigdble=  4.9315e-01 
npeak=  2.5359e-02  area=  l.OOOOe-rOO 

RESULTS 

Effective  support  of  probability  density: 

zsmin=-3 . 5000e-i-00  zsmax=  V.OOOOe+OO  kasmin=  0.0000e+00  kasraax=  3.3C00e-00 

Integrated  value  of  X=  1.0043e+00 

Expectation  value  of  X=  1.0000e+00 

Macroscopic  value  of  X=  2.5359e-02 

Value  X  per  unit  area  =  2.5359e-02 


Figure  2.13:  The  calculation  of  expection  values  of  X. 
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Figure  3.1:  Surface  asperity  in  contact  with  a  rigid  flat  (a  section). 

3  Deformation  Mechanics  of  a  Single  Asperity 

We  now  focus  on  the  analysis  of  a  typical  asperity  in  contact  with  a  rigid  flat.  The  asperity  is 
a  body  of  revolution,  symmetric  about  its  z  =  x3-axis,  and  subjected  to  adhesion  pressures  q 
on  its  exterior  surfaces  that  are  not  in  contact  with  the  rigid  flat,  and  to  contact  pressures  due 
to  its  indentation  into  the  rigid  flat  (see  Fig.  3.1).  The  equations  governing  the  deformation 
of  the  asperity  are  discussed  below. 
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3.1  Momentum  and  Geometric  Equations 

The  momentum  equations  for  the  asperity  are: 


Vijj  =  0 


(3-9) 


where  is  the  Cauchy  stress  tensor  at  a  point  x  =  {x\,xi,x$)  £  (1,  fl  being  the  open 
material  domain  of  the  asperity  and  is  the  divergence  of  the  stress  atj. 

In  rate-dependent  viscoplastic  applications  a  rate  form  of  the  equilibrium  equations  is 
used: 

=  0  (3-10) 

where  the  dot  denotes  the  time  derivative. 

Geometric  equations  express  strains  in  terms  of  displacements: 

1  / 

=  2  (u'u  + 

Strains  can  be  decomposed  into  elastic  and  nonelastic  strains: 


...  _  -(e)  (n) 

ct  j  cij  I  fcjj 


Similarly  as  for  the  momentum  equations,  the  rate  form  of  geometric  equations  will  also 
be  used: 


£>j  =  ^  («.,j  + 


3.2  Constitutive  Equations 

Surface  asperities  for  typical  engineering  surfaces  consist  of  the  same  material  as  the  bulk 
body  with  possible  contaminations  and  structure  change  from  oxidation  and  surface  finish 
processes.  Therefore,  for  general  surfaces  a  variety  of  material  classes  should  be  considered, 
such  as  elastic,  hypoelastic,  elastoplastic,  etc.  In  this  project  two  major  material  classes  are 
considered,  namely: 

•  linearly  elastic  (isotropic  or  anisotropic)  for  some  metal  surfaces,  ceramics,  composites, 
and  hard  rubbers,  and 
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•  viscoelastoplastic  models  for  metallic  surfaces  and  modern  ductile  ceramics. 

3.2.1  Linearly  Elastic  Constitutive  Models 

The  general  linearly  elastic  constitutive  relations  are  given  as: 

<Tij  =  EijklSkl 

where  Exjki  are  the  components  of  the  fourth  order  elasticity  tensor.  It  has  up  to  36  inde¬ 
pendent  coefficients  for  general  anisotropic  materials.  However,  for  most  material  classes, 
the  number  of  material  coefficients  is  much  smaller  and,  for  isotropic  materials,  there  are 
only  two  coefficients,  E  and  u.  The  specific  forms  of  tensor  E  for  various  materials  are  well 
known  and  will  not  be  presented  here. 

3.2.2  Elasto-Viscoplastic  Constitutive  Model  With  Damage 

We  now  describe  the  Bodner-Partom  constitutive  equations  [10,11]  used  in  the  modeling  of 
viscoelastoplastic  asperities.  The  elastic-viscoplastic  analysis  is  based  on  decomposition  of 
strain  rates 


~  +  4n>  (3-11) 

where  superscripts  (e)  and  (n)  denote  elastic  and  nonelastic  strain  components,  respectively. 
The  constitutive  relations  are 


<7  (3.12) 

A  nonelastic  deformation  is  governed  by  the  flow  rule: 

=  fij{<7ij,Zk,Uk) 

Zi  =  ffiio'ij,  Zk) 

u>i  =  hx{<Ttj,u>t) 

where  fij,gx  and  hi  are  constitutive  functions,  z,  are  internal  state  variables,  and  are 
damage  variables.  These  functions  and  state  variables  characterize  the  viscoplastic  response 
of  the  material  with  continuum  damage  effects. 

In  the  particular  version  of  the  Bodner-Partom  theory  applied  in  this  work,  the  nonelastic 
flow  rule  is  of  the  form: 
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4;n)  =  As,, 


u 


where  s,j  are  the  deviatoric  components  of  a  stress  tensor 


"*« j  ‘-b'j  ^  ^kk&ij 


The  current  value  of  parameter  A  is  given  by 


A2  =  Dq  exp 
J  2 


'  z2(l  —  UJ2)\n 


3  J7 


,  A  >  0 


where  J2  is  the  second  invariant  of  a  deviatoric  stress  tensor 


1 


Jt.  — 

Do  is  a  limiting  strain  rate  in  shear,  n  is  a  material  constant  and  2  and  u  are  state  variables 
which  evolve  during  deformation.  In  particular,  2  is  the  hardness  variable,  which  represents 
viscoplastic  hardening  (or  softening)  of  a  material.  The  variable  u  is  the  damage  variable. 
This  variable  represents  weakening  of  the  material  due  to  nucleation  and  propagation  of 
microscopic  voids  and  cracks  in  the  material.  The  micro-cracks  considered  here  are  in  the 
range  of  0.01  mm  in  length.  The  rupture  criterion  is  u;  =  1,  which  corresponds  to  the 
saturation  of  the  material  with  voids.  Alternatively,  a  single  crack  may  grow  to  a  size  on  the 
order  of  1  mm.  In  the  latter  case,  crack  is  too  big  to  be  treated  in  a  continuum  sense,  and 
its  propagation  should  be  followed  using  the  methods  of  fracture  mechanics. 

In  the  framework  of  materials  science  the  value  of  u>  is  usually  interpreted  as  ratio  of  the 
area  of  voids  to  the  total  area  of  a  certain  cross  section  of  a  sample: 


ui  = 


Avoid 


The  state  variables  2  and  u>  evolve  according  to  the  specific  equations  of  the  viscoplastic 
theory: 

1.  Evolution  equations  of  hardness  variable 

The  internal  state  variable  z  consists  of  isotropic  and  directional  components, 


2  =  zl  +  2° 
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The  evolution  equation  proposed  for  the  isotropic  hardening  component  [10,11,24,25]  is 


i;(f)  =  m\[z\ 


zl(t)]Wp(t)-Alzl 


Z!{t)  -  z2 
Z\ 


(3.13) 


with  the  initial  condition,  27(0)  =  z0.  In  the  first  term,  zx  is  the  limiting  (saturation)  value 
of  z1 ,  mj  is  the  hardening  rate,  and  the  plastic  work  rate  is 

Wp  =  <nj :4n) 

which  is  taken  as  the  measure  of  hardening.  22  is  the  minimum  value  of  z 1  at  a  given 
temperature,  and  A\  and  ri  are  temperature  dependent  material  constants.  The  evolution 
form  of  the  directional  hardening  component  (Refs.  [10,11,24,25])  is  defined  as 

Z°(t)  =  fij(t)  Uij(t) 


where  u,j  are  the  direction  cosines  of  the  current  stress  state, 


(3.14) 


The  evolution  equation  for  ftij(t)  has  the  same  general  form  as  that  for  isotropic  hardening 
but  has  tensorial  character, 


where 


and 


Pij  =  m2[zzuij(t )  -  0ij{t)]Wp(t) 


—  A?z\ 


Zx 


v.-j(0 


Vij(t)  =  Pij(t)/[0kl{t)Pkl(t )]» 


Ai(0)  -  0 


As  in  Eq.  (3.13),  m2  is  the  hardening  rate.  A2  and  r2  are  temperature  dependent  material 
constants. 


2.  Evolution  of  damage 
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The  damage  parameter  consists,  in  general,  of  isotropic  and  directional  components, 


u  =  u>!  -f  uD 


The  evolution  of  isotropic  damage  proposed  in  reference  [11]  is  of  the  form 


u1  = 


(3.15) 


In  the  above  P  and  H  are  material  constants,  Q  is  the  stress  intensity  function,  given  by 


Q  = 


+  Bsfih  +  Cl* 


where  is  the  maximum  principal  tensile  stress,  If  is  the  first  stress  invariant  (nonneg¬ 
ative)  and  J?  is  the  previously  introduced  second  invariant  of  deviatoric  stress. 

A,B,C ,  and  u  are  material  constants.  A,B,C  must  satisfy  the  condition 


A  +  B  +  C  =  l 


Clearly,  the  actual  proportion  of  these  constants  selects  the  factor  for  stress  state  which  is 
most  important  in  the  development  of  internal  damage. 

The  initial  condition  for  isotropic  damage  is  u/( 0)  =  0.  In  practical  analyses  the  coef¬ 
ficient  v  is  of  the  order  10  (compare  ref.  [11]).  Thus,  when  SI  (metric)  units  are  used  in 
the  analysis,  the  factor  Q  as  well  as  the  constant  H  reach  extremely  high  values,  beyond  the 
limit  of  real  number  capacity  on  some  computers.  Thus,  for  numerical  analysis,  equation 
(3.15)  was  recast  in  the  equivalent,  but  more  convenient  form: 


where: 

Q  =  Q-  =  Aa+n  +  +  Clf 


P  = 
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The  additional  advantage  of  this  formulation  is  that  both  Q  and  H  are  in  the  stress  units 
(MPa)  instead  of  the  somewhat  cumbersome  ( MPa)u .  The  directional  damage  is  defined  in 
a  manner  very  similar  to  directional  hardening,  namely 


U°  =  uP-Uij 


where  ujj  are  directional  cosines  defined  in  equation  (3.14)  and  the  components  of  a  tensor 
u>D  evolve  according  to  equation 


where  q  and  M  are  material  constants.  The  initial  condition  is 


w?(0)  =  0 


Note  that  there  are  several  problems  with  practical  application  of  directional  damage,  relia¬ 
bility  of  the  above  model  and  conducting  experiments  relevant  for  the  evolution  of  necessary 
parameters.  Even  the  extensive  experiments  presented  in  references  [11,24,25]  did  not  pro¬ 
vide  all  the  necessary  data  and,  hence,  the  damage  model  is  usually  limited  to  the  isotropic 
damage. 


3.3  Boundary  Conditions 

The  asperity  can  be  viewed  as  a  protuberance  of  a  deformable  half  space  (see  Fig.  3.1).  It  is 
subject  to  boundary  conditions  resulting  from  its  support,  contact  with  the  opposing  surface, 
adhesion  and  sliding  resistance. 

3.3.1  Support  Conditions 

If  the  asperity  is  viewed  as  the  protuberance  on  a  deformable  half  space,  the  support  condi¬ 
tions  are  defined  as  zero  displacements  at  infinity: 

lim  u  =  0 

II 2  II  —  oo 
X3<0 

In  practical  computations  we  will  usually  consider  only  a  certain  section  of  the  bulk 
material  surrounding  the  asperity.  Then  the  support  condition  will  be: 

u  =  0  on  r„ 


on  the  cut-off  boundary  fu. 


45 


3.3.2  Contact  Condition 


Let  the  position  of  the  rigid  flat  (see  figure  3.2)  be  defined  by: 

•  a  point  p0(x0,y0,  z0)  which  belongs  to  the  flat,  where  x0,y0,z0  are  its  coordinates  in 
the  initial  configuration, 

•  unit  vector  N,  normal  to  the  flat, 

•  displacement  w  of  the  flat  in  direction  N. 

Separation  of  material  point  in  the  deformable  body  from  the  flat  is  then  given  by  the 
following  formula 


s  =  (x  +  u  -  p0)  ■  N  —  w  —  d 

where: 

x  -  initial  position  of  a  material  point, 
u  -  displacement  of  this  point, 
d  -  intermolecular  distance  which  is 
important  when  adhesion  is  taken 
into  account,  otherwise  d  =  0. 

The  condition  that  the  asperity  cannot  penetrate  the  rigid  flat  is: 

s  >  0  on  T 

The  actual  contact  region  is  Tc  =  {x  £  T,  s(a:)  =  0}.  The  difficulty  associated  with  the 
contact  condition  in  the  above  form  is  that  it  results  in  a  weak  formulation  of  the  problem 
in  the  form  of  variational  inequality,  rather  than  the  equation.  In  order  to  avoid  difficul¬ 
ties  involved  in  solving  variational  inequalities,  the  contact  condition  is  usually  regularized 
[54,71].  In  this  work  we  will  use  the  penalty-type  regularization  of  the  form: 

=  ^v(°) on  rc 

where  a  =  —  s  is  the  approach  (penetration)  and  tcN  is  the  value  of  traction  normal  to 
the  flat  which  defines  resistance  of  the  surface  to  penetration.  Because  we  will  be  using 
rate  formulation  in  viscoelastoplastic  analysis,  it  is  beneficial  to  introduce  a  continuously 
differentiable  penalty  function,  for  example  in  the  form  presented  graphically  in  Fig.  3.2: 
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Figure  3.2:  Penalty  function  for  contact  condition. 
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Here  H  is  a  large  number  (normal  stress)  and  £  is  a  small  number  (penetration).  The 
above  penalty  function  guarantees  continuous  derivative  of  the  normal  traction  with  respect 
to  a,  which  greatly  improves  practical  performance  of  the  numerical  computations. 
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3.3.3  Adhesion 


An  important  contributing  factor  to  friction  on  contact  surfaces  is  adhesion:  the  intermolec- 
ular  attractive  forces  that  depend  on  atomic  spacing  and  the  corresponding  surface  energies 
of  materials.  For  highly  polished  uncontaminated  surfaces,  adhesion  forces  can  be  very  large, 
leading  to  the  virtual  welding  of  one  surface  to  another,  while  for  rough  contaminated  sur¬ 
faces,  adhesion  effects  are  often  negligible.  For  engineering  surfaces  under  common  working 
conditions,  adhesion  effects  can  be  significant,  so  that  a  rational  model  of  contact  and  friction 
should  take  them  into  account. 

The  DMT  adhesion  model,  proposed  by  Derjaguin,  Muller,  and  Toporov  [36]  and  later 
refined  by  Muller,  Derjaguin,  and  Toporov  [66],  attempts  to  characterize  the  attractive  forces 
on  a  sperical  elastic  asperity  in  contact  with  a  rigid  flat,  assuming  that  the  shape  of  the 
deformed  asperity  is  given  by  the  Hertz  theory  and  that  no  attractive  forces  exist  in  the 
contact  region.  The  JKR  model,  due  to  Johnson,  Kendall,  and  Roberts  [52]  also  analyzes 
the  elastic  spherical  asperity-fiat  problem  with  Hertz  theory,  but  assumes  that  attractive 
forces  are  confined  to  the  contact  area  and  that  the  attractive  forces  produce  an  elastic 
deformation  of  the  asperity.  The  JKR  model  has  been  found  to  be  more  suitable  for  soft 
materials,  such  as  rubbers,  while  the  DMT  model  is  claimed  to  be  more  suitable  for  harder 
materials  with  high  surface  energies  (see  Chang,  Etsion.  and  Bogy  [26-28]  and  Pashley 
[73]).  Survey  papers  on  developments  in  adhesion  models  have  been  contributed  by  Pashley, 
Pethica,  and  Tabor  [75]  and  bv  Pashlev  and  Pethica  [74].  See  also  MacFarlane  and  Tabor 

[58]. 

We  shall  include  adhesion  effects  in  our  contact  and  friction  theory  by  using  an  approach 
similar  to,  but  more  general  than  that  of  the  DMT  model.  We  continue  to  assume  that  the 
surfaces  are  isotropic,  rough,  and  have  a  Gaussian  distribution  of  peak  heights,  that  there  is 
no  interaction  between  asperities,  and  that  it  suffices  to  consider  a  single  asperity  impending 
on  a  rigid  flat.  Following  Muller,  et  al.  [67],  we  characterize  the  attractive  adhesion  pressure 
q(s)  as  that  attractive  force  per  unit  surface  area,  acting  normal  to  the  mean  asperity  height 
plane,  resulting  from  the  Lennard-Jones  interaction  (interatomic)  potential  <j>  of  surface 
physical  chemistry.  Then, 


where 


(3.16) 


s  =  the  separation  of  the  two  surfaces  outside  the  contact  area 
d  =  the  intermolecular  distance  (generally  d  =  0.3  —  0.5nm) 
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In  (3.6),  A7  is  the  surface  energy  of  adhesion  and  is  defined  as  follows:  if  71  and  72  are 
the  surface  energies  of  surfaces  1  and  2,  respectively,  before  contact,  and  712  is  the  joint 
surface  energy  of  the  interface  after  contact,  then 


A7  =  71  +  72  -  712  (3.17) 

For  values  of  surface  energies  of  adhesion  for  various  metals,  see  Rabinowicz  [78]  or 
Ferrante,  Smith,  and  Rose  [39].  Traction  resulting  from  adhesion  can  be  expressed  by  the 
following  formula 


ta  =  -q{s)N  =  q{a)N 

where  s  is  calculated  according  to  the  formula  defined  in  the  previous  subsection. 

The  principal  mathematical  difficulty  inherent  in  charactrizing  the  adhesion  pressure  is 
that  it  is  developed  only  on  surface  material  outside  the  contact  area,  which,  a  priori ,  is 
unknown.  Several  concluding  remarks  on  adhesion,  however,  are  in  order  at  this  point. 


Remarks 


1. 


Fuller  and  Tabor  [42],  using  the  JKR  model  of  adhesion  and  the  rough  surface  asperity- 
based  model  of  Greenwood  and  Williamson  [44],  presented  an  adhesion  parameter  9  of 
the  form 


A7  V  R 

where  E  is  the  effective  modulus  of  elasticity  of  the  contact  surfaces 


(3.18) 


E-'  =  (l  -  „J)  Ep  +  (l  -  vl)  E~P 


Ej,V{  being  Young’s  modulus  and  Poisson’s  ratio  of  surface  i,  a  is  the  standard  de¬ 
viation  of  surface  heights,  and  R  the  mean  radius  of  spherical  elastic  asperities.  The 
larger  the  value  of  9 ,  the  less  significant  the  effects  of  adhesion,  and  for  rubber  spheres, 
experiments  showed  that  adhesion  become  negligible  for  9  >  10. 


2.  Chang,  Etsion,  and  Bogy  [26-28]  presented  a  study  of  adhesion  effects  using  the  DMT 
model  and  an  elasto-plastic  asperity  model  of  the  contact  surface.  They  investigated 
the  importance  of  adhesion  with  varying  values  of  surface  energy  and  plasticity,  as 
measured  by  the  plasticity  index  of  Greenwood  and  Williamson  [44], 


2  E  fa 
ttKH  V  R 


(3.19) 
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Figure  3.3:  Pressure  on  a  deformed  asperity;  the  adhesion  pressures  q  due  to  molecular 
attraction  and  the  pressure  p  on  the  contact  surface. 
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where  H  is  the  Brinell  hardness  (of  the  softer  of  the  two  materials)  and  K  =  0.454  + 
0.41j/i,i/i  being  the  Poisson  ratio  of  this  material.  These  authors  conluded  that  the 
“pull-off  force”  due  to  adhesion  (i.e.,  the  integrated  suction  force  due  to  adhesion) 
becomes  negligible  for  hard  steel  when  the  adhesion  index  9  is  greater  than  100,  as 
compared  to  9  >  10  for  rubber;  the  adhesion  force  is  negligible  compared  to  the  contact 
load  when  the  plasticity  index  $  >  2.5  or  when  the  surface  energy  A  7  <  0.5  ;/m2 
for  a  sufficiently  small  external  force.  They  concluded  that  “for  smooth  clean  surfaces 
the  adhesion  can  be  well  over  20  percent  of  the  contact  load  and,  thus,  cannot  be 
neglected.” 

3.  Adhesion  forces  are  time  dependent  and  generally  increase  with  time  of  contact,  even¬ 
tually  acquiring  a  constant  value  for  static  contact.  Thus,  the  static  adhesion  models 
generally  attempt  to  predict  the  maximum  value  of  adhesion  forces  and  can  overesti¬ 
mate  adhesion  effects  for  dynamic  contact. 

3.3.4  Shear  Resistance 

To  construct  new  constitutive  models  of  friction  it  is  necessary  to  characterize  the  resistance 
of  the  rough  interface  to  sliding  (i.e.,  to  tangential  motions  of  the  reference  planes  relative 
to  one  another).  While  this  aspect  of  the  modeling  approach  still  requires  much  study,  there 
appears  to  be  at  least  three  methods  available  for  this  purpose.  First,  Bowden  and  Tabor 
[13,  14]  estimated  the  resistance  to  impending  motion  (more  precisely,  the  static  coefficient 
of  friction)  by  calculating  the  shear  strength  of  metallic  oxide  junctions  developed  on  the 
contact  surface.  Similarly,  Chang,  Etsion,  and  Bogy  [26-28]  calculated  the  tangential  load 
required  to  reach  the  fracture  strength  of  metallic  junctions  as  an  indication  of  the  tangential 
force  required  to  produce  sliding.  Villiaggio  [97],  on  the  other  hand,  studied  the  problem  of 
contact  of  periodically  spaced  elastic  asperities  and  defined  the  load  at  which  sliding  initiates 
as  that  which  reduces  the  curvature  of  the  resisting  elastic  asperities  to  zero.  Francis  [41] 
modeled  the  micro  sliding  resistance  using  empirical  relations  based  on  existing  experimental 
data. 

Extensive  sequences  of  experiments  on  sliding  resistance  of  thin  films,  involving  27  dif¬ 
ferent  materials  and  a  wide  range  of  normal  loads,  are  described  in  the  papers  of  Boyd  and 
Robertson  [15],  Briscoe,  Scruton,  and  Willis  [17],  Towle  [92],  and  Briscoe  and  Tabor  [18].  In 
all  of  these  studies,  it  was  discovered  that  (on  a  microcontact  interface)  the  interfacial  shear 
stress  t3  during  sliding  was  a  function  of  the  normal  stress  crn  =  crun;n,,  and,  according  to 
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Francis  [41],  their  empirical  findings  suggest  that 


for  light  loads 

for  intermediate  loads 

for  heavy  loads 


T, 


ci  +  c2<rn 

c3<t™  (0.6  <  m  <  1.4) 

'  log-1  [c5  +  c(an)log<7„] 

where  0.0  <  c(crn)  <  1.9  ;  ^  >  0 

d(logcrn) 


(3.20) 


wherein  loads  were  varied  over  a  factor  of  10  or  more  and  starting  from  15  MPa  (light),  40 
MPa  (intermedite)  and  200  MPa  (heavy).  Francis  [41]  points  out  that  a  good  approximation 
to  all  of  these  cases  is  the  simple  quadratic  function, 


Ts  —  Cq  +  Ci<7„  +  c2cr^ 


(3.21) 


where  Co,C!,c2  are  material  constants. 

Once  a  micro-shear  resistance  is  characterized,  the  macro  sliding  resistance  can  be  com¬ 
puted  using  the  statistical  summation  procedures  described  earlier. 

3.3.5  Initial  Conditions 

Smooth  functions  u°(a;)  and  z°{x)  are  prescribed  such  that  for  x  €  Q, 


t*<(*i0)  =  «?(*) 

Zi{x,  0)  =  2°(x) 
Wi(sc,0)  =  0 
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Figure  3.4:  A  regularization  function  <f> 

3.4  Variational  Formulation 

In  order  to  obtain  a  weak  formulation  of  a  boundary-value  problem  we  introduce  the  space 
of  functions 

V  =  {v  €  ,  v(x)  ->  0  as  |[x||  ->  oo} 

where  Cl  is  a  computational  domain,  N  is  the  dimension  of  the  physical  space  (2  or  3),  and 
VFm’p(fl)  is  the  Sobolev  space,  where  specific  values  of  m,p  and  q  depend  on  the  particular 
form  of  constitutive  equations. 

Multiplying  the  equilibrium  equation  (3.10)  by  a  test  function  and  integrating  over  Cl  we 
obtain  the  weak  form  of  the  rate  equilibrium  equations: 

ViVijj  dCl  =  0  V  v  £  V 

After  the  substitution  of  the  constitutive  equations,  application  of  the  divergence  theorem 
and  a  grouping  of  terms  the  following  variational  problem  is  obtained: 
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Find  a  displacement  rate  field  t  — >  u(x,  t)  €  V  such  that 

f  t >ijEijUukjdn  =  [  VijEijuitfdn  +  /  U.t.ds  V»  €  V  (3.22) 

jn  -'an 

The  rates  of  nonelastic  strains  can  be  obtained  from  the  relevant  constitutive  theory 
(see  Section  3.2.2).  For  elastic  materials  they  are  identically  equal  to  zero.  Therefore  for  these 
materials  it  is  also  possible  and  computationally  more  efficient  to  use  a  total  formulation  of 
the  problem  which  is  as  follows. 

Find  a  displacement  field  u(x)eV,  such  that 

/  Vi  ,Eijkiiik  idfi,  =  f  Vit{ds  V  v  €  V  (3.23) 

Ja  '  Jd n 

Note  that  the  values  or  the  rates  of  tractions  on  <90  need  to  be  expressed  in  terms  of 
displacements  using  formulas  presented  in  previous  sections  (contact  condition,  adhesion, 
and  sliding  resistance). 

3.4.1  Boundary  Integrals 

The  boundary  <90  can  be  decomposed  in  the  following  way 

<90  =  ru  u  r£  u  ra  u  rc 


where 


ru  —  support  zone  (kinematic  boundary  conditions), 

T,  —  static  zone  (static  boundary  conditions), 

Ta  —  adhesion  zone  (adhesion  traction  is  not  negligible), 

Tc  —  contact  zone(a  >  0). 

The  integrals  over  dfl  in  formulas  (3.22)  and  (3.23)  can  be  calculated  as  sums  of  integrals 
over  these  four  parts  of  the  boundary. 

*  Support  zone 

According  to  the  formula  (3.3.1) 


itj  =  0  on  ru 
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rigid  flat 


Figure  3.5:  Separation  and  contact  when  intermolecular  forces  are  taken  into  account. 

so  there  is  no  need  to  compute  integrals  over  Tu 

•  Static  zone 

Traction  on  this  part  of  the  boundary  is  known,  so  it  can  be  integrated.  In  the  case  of 
asperities  this  traction  is  usually  equal  to  zero. 

•  Adhesion  zone 

The  integrals  over  Ta  have  the  following  forms: 

(a)  total  formulation 


(b)  rate  formulation 
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where 


ViawNids 


f-V 

\a  ) 

When  adhesion  forces  are  taken  into  account  then  the  intermolecular  distance  d  plays  an 
important  role.  A  graph  of  traction  resulting  from  adhesion  for  d  =  0.5  nm  and  A7=  0.1 
J  /  m2  in  terms  of  the  separation  is  shown  in  figure  3.6.  Clearly,  the  adhesion  traction  is 
strongly  nonlinear,  and  practically  vanishes  at  separtions  greater  than  4d. 

Integration  of  adhesion  derivative  with  respect  to  separation  is  required  in  formulas  3.4.1. 
Figure  3.7  shows  this  function.  Because  the  function  changes  very  rapidly,  it  has  to  be 
integrated  with  higher  accuracy  than  e.g.  shape  functions.  Some  numerical  experiments 
were  carried  out  to  check  how  many  integration  points  are  necessary  to  compute  accurately 
these  integrals.  They  showed  that  10  integration  points  in  one  direction  are  enough  to 
integrate  influence  of  adhesion  forces  over  an  element  face  which  has  dimensions  less  than 
10*d,  providing  that  the  slope  i.e.  angle  between  the  face  and  the  flat  is  small  (less  than 
10°).  We  have  such  a  situation  in  the  case  of  asperities.  Moreover,  integration  need  only  be 
performed  on  part  of  ro  on  which  separation  is  less  than  A  *  d.  At  this  and  higher  distance, 
adhesion  traction  and  its  derivatives  are  practically  equal  to  zero. 

•  Tractions  on  the  contact  zone  rc 

Penalty  method  which  we  use  here  to  solve  variational  inequality  problem  is  equivalent 
to  allowing  for  penetration  of  the  rigid  flat  by  the  asperity  (a  >  0).  But  this  penetration 
results  in  a  big  normal  traction  in  the  form 

1 

t/vr  —  — aiV 

t 

where  e  is  a  small  number. 

In  addition  to  normal  traction  there  exists  also  a  tangential  traction  (sliding  resistance) 
on  the  contact  zone.  There  exist  several  theories  which  give  formulas  for  sliding  resistance 
in  terms  of  normal  pressure  [41]  see  section  4.3.4  for  discussion.  Generally  all  of  them  can 
be  expressed  in  the  following  form 


t  = 


0  when  6  =  0 

when  6^0 


(3-24) 
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derivative  of  adhesion  [MPa/m] 


X  Graph 


where 


t  —  sliding  resistance, 

h  —  function  which  defines  the  relation  between 

sliding  and  normal  tractions,  its  form  depends 
on  adopted  model  of  sliding  resistance, 

6  —  vector  of  relative  displacement  between 
the  asperity  and  the  flat. 

Components  of  t  are  not  continuous  with  respect  to  6  when  6  =  0.  To  facilitate  numerical 
solution  of  the  contact  problem  a  regularization  of  the  function  t  is  introduced.  The  regu¬ 
larization  will  be  done  for  components  of  sliding  resistance.  This  vector  has  at  most  two  non 
zero  components.  Let  the  first  component  t$  have  direction  of  an  arbitrary  vector  5  parallel 
to  the  flat.  Let  the  second  component  tj  have  direction  of  vector  product  iV  x  S  =  T,  where 
N  is  normal  to  the  flat  and  pointing  towards  the  asperity.  If  6s  and  br  are  coordinates  of  6 
then  the  components  of  t  after  regularization  can  be  expressed  by  the  following  formulas 

ts  =  6(t^)^(6)^s(6s,  6j) 
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tr  =  M*N)<£(i>)'J,r(&s>6r) 


where  <f>  is  a  function  which  provides  smooth  approximation  (of  class  C1)  for  function  h, 
and  has  the  following  form 


-1 

for 

*  <  -4 

2— 

<6 

+  (*)’ 

for 

— 4  <  x  <  0 

2— 

-te) 

for 

0<x<eb 

1 

for 

eb  <  x 

(3.25) 


where  eb  is  a  small  number. 

tfs  and  tyj  are  functions  which  guarantee  a  proper  decomposition  of  vector  £,  i.e.  they 
provide  that  ||t||  =  t2s  +  t\.  Functions  have  the  following  forms 


#s  = 


bs 

\jb2s  +  b\ 


— 


bf 


Combination  of  the  above  formulas  leads  to  the  following  contributions  to  the  weak 
statement  of  the  problems: 

for  the  total  formulation 


f  Vitids  =  f  Vi(-Ni  +  tsSi  +  tTTi)ds 
J rc  J rc  c 

for  the  incremental  formulation 

J  v^ids  =  J  ViCtijitjds  + 

+ £  v,  [u.v, + ^ * + (If" + r‘ 

where 


ds 


Q'i 


-  ;W,  +  ^4£w,  +  f£sir,+ 


(3.26) 
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(3.27) 


dtr 

8a 


TiN‘  +  m?-s>  + 


dtx 

dbr 


T{Tj 


w  —  displacement  of  the  flat  in  N  direction, 

6  —  displacement  of  the  fiat  in  S  direction. 

3.5  Solution  Method  for  Elastic  Contact  Problems 

Formulation  of  the  contact  problem  is  nonlinear  even  in  the  case  of  contact  with  an  elastic 
body  because  the  area  of  contact  depends  on  displacements  (rc  =  Tc(w)).  Generally,  requir¬ 
ing  that  variational  equation  (3.23)  be  satisfied  for  every  test  function  leads  to  the  following 
nonlinear  system  of  equations: 


L(u)~R  =  0,  (3.28) 

where  L  stands  for  the  left-hand  side  and  R  for  the  right-hand  side. 

To  solve  the  problem  effectively,  Newton- Raphson  iteration  technique  was  used.  The  idea  of 
the  method  is  to  substitute  a  nonlinear  functional  by  its  linear  part.  Linearization  is  made  at 
a  series  of  points.  Each  point  is  the  solution  of  the  problem  linearized  at  the  previous  point. 
If  the  series  is  convergent,  it  is  convergent  to  a  solution  of  the  nonlinear  problem.  For  a 
simple  case  of  a  simple  nonlinear  equation  with  one  unknown,  two  steps  of  Newton-Raphson 
method  are  shown  in  Fig.  3.8. 

Basic  formulas  of  the  Newton-Raphson  method  are  presented  below.  Let  us  assume  that 
we  know  a  field  un  which  is  an  approximate  solution  of  the  equation  (3.28).  u0  =  0  can  be 
assumed.  First  two  components  of  the  Taylor  series  evolution  of  left  hand  side  of  equation 
(3.28)  give 


L(un)  —  R  —  gradw[L(un)  —  R]6u  ss  0 
Assuming  that  8u  —  i*n+1  —  un  we  obtain 

L{un)  -  R  +  gradttL(tin)(un+1  -  un)  =  0 

Equation  (3.28)  is  a  linear  equation  for  un+l.  We  use  this  equation  to  compute  a  sequence 
of  approximate  solutions  u ',...,uM.  The  process  is  stopped  when  Hu^  —  uM-1||  and 
||£(uAf)  —  i2||  are  small  enough. 

For  contact  with  elastic  bodies,  equation  (3.28)  has  the  following  form 
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J^VijEijklunkjldn+Jr  vtatJu"+lds+ VifajU1-*1  ds  =  u,[a,-,ujn)  +pi}ds+ u, [/?,_, ujn)  +q,]ds 

where 


C*ij 

Pa  = 

Pi  = 

<7.  = 

Variations  were  evaluated  neglecting  dependence  rc  =  rc(rt).  The  above  linearized  prob¬ 
lem  is  solved  by  the  standard  FEM. 

In  order  to  provide  automatic  control  of  the  performance  of  nonlinear  procedures,  an 
expert  system-like  approach  has  been  applied.  This  application  is  based  on  our  previous 
research  on  automation  of  computational  procedures  [96],  and  employs  several  heuristic 
rules  to  monitor  and  control  the  performance  of  nonlinear  iterations.  While  in  the  original 
implementation  discussed  in  reference  [96]  the  specialized  knowledge  engineering  software 
was  used  to  develop  the  expert  system,  in  this  project  the  essential  features  of  the  expert 
system  were  coded  in  FORTRAN  and  included  in  the  code.  The  expert  system  is  activated 
at  each  time  step  after  completing  a  prescribed  number  of  iterations  (sufficient  to  estimate 
trends  in  error  histories).  The  decisions  of  the  expert  system  are  used  to  control  the  solution 
process  and  obtain  a  converged  solution  at  minimum  cost. 

3.6  Solution  Method  for  Viscoplastic  Contact  Problems 

Formulation  of  the  problem  in  this  case  is  time-dependent. 

The  strategy  employed  in  the  solution  of  this  problem  is  as  follows:  with  the  initial 
distribution  of  stress,  temperature  and  internal  variables  specified  use  the  rate  form  of  the 
equilibrium  condition  (Eq.  (2.23))  to  obtain  the  nodal  displacement  rates.  Then  integrate 
the  constitutive  equations  forward  in  time  at  the  element  Gauss  integration  points.  With 
updated  value  of  the  stress,  temperature  and  internal  variables  at  the  new  time,  the  equilib¬ 
rium  equation  is  solved  again.  This  sequence  of  determining  the  nodal  displacement  rates, 


N,Nj 


is  the  same  as  in  the  incremental  formulation 

sM3(fT-  "'31 

da  \a ) 

-(w  +  d-  PN  +  P0N)N, 

3  /d'91 


8  A7 

3  d 


'd 

.  a , 


,aJ 


Ni 


(3.29) 

(3.30) 

(3.31) 

(3.32) 
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then  advancing  the  constitutive  equations  in  time  is  continued  until  the  desired  history  of 
the  initial  boundary-value  problem  has  been  obtained. 

Thus,  the  algorithm  proceeds  through  the  following  steps: 

1.  At  time  f,  initialize  cr^-,  Z;  for  each  element; 

2.  Calculate  Zk)  for  each  element; 

3.  Assemble  and  solve  [K]U  =  F\ 

4.  Calculate  for  each  element,  e  =  [B\U\ 

5.  Calculate  cr,j  for  each  element,  a  —  [£](e  —  en); 

6.  Calculate  Z,  for  each  element,  Z,  =  Zk ); 

7.  Integrate  djj,  Z  forward  for  each  element  to  get  <7,j  and  Z,  at  i  +  A t3; 

8.  If  t  -f  At,  <  go  to  2,  otherwise  stop. 

The  computational  method  above  has  been  presented  for  a  constant  time  step  A <s- 
Computational  experience  by  several  investigators  (see  refs.  [7,8,55])  indicates  that  a  very 
small  time  step  can  be  required  because  of  the  “stiff’  nature  of  the  ordinary  differential 
equations  describing  the  internal  state  variables.  To  gain  improved  efficiency  and  reliability 
a  variable  time  step  algorithm  has  been  implemented.  The  basic  idea  of  this  variable  time 
step  algorithm  is  presented  below  for  a  scalar  evolution  equation. 

The  solution  is  advanced  using  a  predictor- corrector  scheme.  The  predictor  phase  consists 
of  an  Euler  step: 


■4^ 

II 

(3.33) 

The  solution  is  advanced  using  a  predictor-corrector  scheme.  The  predictor  phase  consists 
of  an  Euler  step: 

yf+a  t  =  vt  +  Myt 

(3.34) 

ih+At  ~  f  {yt+Ati*  "h 

An  error  indicator  E  [15,24]  is  then  computed  from 

(3.35) 
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(3.36) 


p  ..  (yS-At-yt)l 

2|l/j+Ai| 

The  error  indicator  is  next  compared  with  a  preset  error  criterion  and  if  the  criterion  is 
met,  the  time  step  is  small  enough  to  proceed  to  the  corrector  stage.  Otherwise,  the  predictor 
phase  for  Eqs.  (3.34)-(3.35)  is  repeated  with  a  smaller  time  step.  For  the  viscoplastic 
evolution  equations  with  damage  modeling,  the  control  variables  used  to  calculate  the  error 
indicator  were  the  components  of  a  stress  tensor  cr,j,  internal  state  variables  Z,,  and  the 
damage  variables  u>,-,  with  the  maximum  of  these  selected  as  the  controlling  error. 

The  corrector  phase  is  the  modified  Newton  scheme, 

llavg  =  (j/t  +  2/j+At)  /2 
Vt+At  ■“  Vt  "b  Ati/avg 

A  flowchart  depicting  the  adaptive  scheme  is  shown  in  Fig.  3.9. 
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Figure  3.9:  Flowchart  for  the  solution  of  viscoelastoplastic  evolution  problem  with  adaptive 
timestepping 
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4  Finite  Element  Analysis  of  Contact  Problems  with 
Friction 

4.1  General  Information  About  the  3D  Finite  Element  Code 

The  development  of  asperity  modeling  capabilities  was  based  on  existing,  state-of-the-art 
three-dimensional  adaptive  finite  element  kernel  code,  which  consists  of  several  separate 
modules  organized  around  the  common  data  structure  and  execution  supervisor.  Figure  4.1 
shows  a  general  structure  of  the  code  directories.  The  most  important  modules  are: 

•  an  object-based  data  structure  designed  specifically  for  the  h-p  adaptive  finite  element 
method, 

•  an  execution  supervisor  controlling  the  overall  execution  of  the  computations, 

•  pre-  and  postprocessors, 

•  an  adaptive  package, 

•  linear  equation  solvers,  and 

•  a  solver  for  a  specific  boundary  value  problem  (in  this  case,  asperity  modeling). 

Importantly,  the  elements  of  the  finite  element  data  structure,  adaptive  package,  graphical 
interfaces  and  linear  equations  solvers  were  designed  to  be  applicable  to  a  general  class  of 
problems,  and  relatively  easy  customizable  to  specific  problems  in  solid  mechanics  or  fluid 
mechanics.  Below,  selected  modules  of  the  above  kernel  are  discussed  in  more  detail. 

Object-based  Data  Structure 

A  new  state-of-the-art  data  structure  was  designed  and  implemented  in  the  kernel  to  avoid 
typical  limitations  of  traditional  finite  element  codes,  such  as: 

•  fixed  size  common  blocks  and  arrays, 

•  predefined  limits  on  problem  size, 

•  element  information  spread  throughout  memory  in  variety  of  arrays. 

The  object-based  data  structure  was  coded  in  C  computer  language,  which  allows  for 
dynamic  memory  allocation  and  more  flexible  handling  of  objects  and  structures.  Typical 
examples  of  objects  handled  by  this  data  structure  are: 
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•  elements, 

•  nodes, 

•  boundary  condition  data, 

•  set  of  degrees  of  freedom,  etc. 

The  major  advantages  of  object-based  handling  of  these  structures  are  listed  below: 

•  the  objects  are  created  only  when  needed, 

•  all  related  information  is  contained  in  one  structure,  and  closely  packed  in  memory, 

•  all  the  objects  are  automatically  saved/restarted, 

•  the  memory  is  reused  when  object  is  deleted,  say  during  mesh  refinement. 

It  is  of  importance  to  note  here,  that  the  elements  in  the  above  data  structure  are 
grouped  and  colored  in  order  to  facilitate  vector  and  parallel  processing.  The  basic  idea  of 
this  vectorization  and  parallelization  is  presented  in  figure  4.2.  Importantly,  all  the  elements 
in  a  batch  are  of  the  same  type,  so  that  the  generation  of  element  stiffness  matrices  and 
right-hand  sides  can  be  effectively  vectorized  by  putting  loop  over  elements  as  the  innermost 
loop.  On  the  other  hand,  since  the  elements  in  different  colors  have  no  common  nodes  or 
sides,  the  generation  of  element  of  element  matrices  and  assembly  for  different  colors  can  be 
performed  in  parallel. 

Adaptive  three-dimensional  finite  element  meshes 

The  finite  element  kernel  is  designed  to  handle  h-p  adaptive  finite  element  meshed  for 
three-  and  two-  dimensional  problems.  By  h-p  adaption  we  understand  a  finite  element 
technique,  wherein  the  elements  can  be  automatically  subdivided  into  smaller  elements  (h- 
refinement  /unrefinement)  and  the  polynomial  order  of  approximation  can  be  locally  in¬ 
creased  or  reduced.  A  major  advantage  of  properly  designed  h-p  mesh  is  that  it  can  achieve 
a  higher  order  of  accuracy  with  much  less  degrees  of  freedom  than  traditional  finite  element 
methods.  Moreover,  the  optimal  mesh  is  designed  automatically  by  adaptive  procedure 
driven  by  appropriate  error  estimators. 

For  three-dimensional  problems,  anisotropic  h- refinement  can  offer  a  wide  improvement  in 
computational  effort  over  more  conventional  isotropic  refinement  schemes.  This  is  primarily 
true  because  anisotropic  refinement  allows  for  selected  refinement  in  the  directions  of  interest 
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batch 


coior  —  not  connected  elements 

group  — "topologically''  identical  elements 

batch  —  optimal  vector  length 


Figure  4.2:  Vectorization  and  parailelizaton  for  groups  and  colors 


only  (i.e.,  directions  of  high  error).  Thus,  anisotropic  refinement  may  greatly  reduce  the  total 
number  of  unknowns  in  many  problems,  in  turn  reducing  the  required  computational  effort. 

Isotropic  refinement  implies  that  an  element  is  identically  refined  in  each  local  direction. 
For  a  hexahedral  element,  an  isotropic  refinement  is  a  division  into  two  along  each  of  the  three 
local  directions,  which  results  in  eight  sub-elements.  In  contrast,  an  anisotropic  refinement 
of  a  hexahedral  element  is  a  division  into  two  along  a  single  local  direction,  resulting,  of 
course,  in  only  two  sub-elements.  Thus,  if  solution  phenomena  is  oriented  with  respect  to 
a  particular  local  direction,  then  anisotropic  refinement  allows  for  degrees  of  freedom  to  be 
introduced  only  in  the  direction  which  actually  reduces  the  total  error.  Isotropic  refinement, 
on  the  other  hand,  would  have  introduced  degrees  of  freedom  in  all  directions,  many  of  them 
providing  little  improvement  to  the  overall  solution.  Anisotropic  refinement  can,  therefore, 
provide  a  higher  level  of  accuracy  than  isotropic  refinement  using  the  same  number  of  degrees 
of  freedom. 

Several  examples  of  h-adapted  meshes  in  three  dimensions  will  be  shown  in  Section  6. 
Note  that  the  mesh  refinement  introduces  several  theoretical  and  numerical  complications 
into  the  algorithm,  such  as: 

•  constrained  or  ’’hanging”  nodes  between  elements  of  different  refinement  level, 

•  propagation  of  constraints  and  possible  ’’deadlocks”  in  the  case  of  directional  refine¬ 
ments  for  complex  geometries, 

•  complications  of  unrefinement  due  to  one-to-two  approximation  rule. 

Detailed  discussion  of  these  issues  is  beyond  the  scope  of  this  report.  It  is  sufficient  to 
note,  that  before  application  of  the  above  kernel  to  asperity  modeling  all  these  difficulties 
have  been  successfully  resolved  and  the  existing  kernel  offers  operational  unique  automated 
directional  refinement  capability  for  three-dimensional  hexagonal  meshes. 

Interactive  user’s  interface  and  graphical  postprocessing 

have  been  implemented  in  the  adaptive  kernel  to  enable  user-friendly  operation  of  the  code 
and  viewing  of  three-dimensional  result.  The  interactive  graphic  interface  is  based  on  a 
window  environment,  with  a  menu-driven  selection  of  options.  A  sample  view  of  the  screen 
with  several  windows  open  is  presented  in  figure  4.3. 

The  graphic  interface  can  be  customized  for  specific  applications,  such  as  contact  and  fric¬ 
tion  modeling,  so  that  the  solution  process  and  essential  data  can  be  controlled  interactively 
by  the  user. 
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Figure  4.3:  Sample  screen  with  interactive  window  environment 
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The  most  important  feature  of  the  graphical  user  interface  is  a  three-dimensional  inter¬ 
active  postprocessing  capability.  For  two-dimensional  models,  such  visualization  is  rather 
trivial  as  all  of  the  computational  domain  is  always  visible  and  it  is  a  simple  matter  to  zoom 
and/or  pan  through  the  mesh  to  closely  review  the  results.  The  real  challenge  comes  from 
the  need  to  visualize  phenomena  in  three-dimensional  domains  where  most  of  the  numerical 
data  is  actually  hidden  from  the  observer  and  one  needs  to  enter  the  domain  to  view  the 
local  structure  of  the  solution. 

The  postprocessing  capability  implemented  in  the  kernel  is  capable  of  displaying  solution 
obtained  on  structured  and  unstructured  meshes,  with  both  h-refinement  and  p-enrichment 
present  in  the  mesh.  The  package  is  fully  interactive  and  operates  efficiently  on  high-end 
workstations.  The  basic  graphic  features  displayed  include: 

•  mesh  plots, 

•  isosurfaces  of  selected  quantities, 

•  slicing  planes  with  overlaying  isolines, 

•  deformed  configurations, 

•  three-dimensional  cursor  for  picking  pointwise  values  of  the  solution, 

All  the  above  displays  are  available  with  interactive  translation,  rotation  and  zoom  op¬ 
tions,  hidden  line  removal,  panning,  etc. 

Importantly,  the  graphics  package  is  designed  to  take  advantage  of  specialized  graphic 
hardware  and  software  available  on  many  platforms.  The  primary  platform  for  the  package 
is  the  SGI  Iris  family,  which  is  also  a  primary  platform  in  this  project.  Alternatively,  X- 
windows  graphics  is  supported,  which  is  operational  on  most  Unix  workstations. 

4.2  Formulation  of  a  Structural  Deformation  Problem  in  the  3D 
Code 

In  the  project,  the  3D  kernel  was  customized  to  solve  contact  problems.  It  was  supplemented 
with  over  8,000  instructions.  They  enable  to  run  specialized  drivers  when  a  contact  problem 
is  to  be  solved.  The  drivers  solve  the  contact  problem  using  either  total  formulation  or 
incremental  formulations  (see  section  3.4).  The  first  one  uses  Newton- Raphson  iterative 
method  and  calls  the  FEM  linear  solver  at  each  iteration  step.  The  second  driver  uses  Euler 
predictor  corrector  integration  method  with  automatic  time  step  control  and  calls  the  FEM 
linear  solver  twice  at  each  time  step. 
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To  solve  solid  mechanics  problems  with  contact,  additional  customization  of  the  kernel 
FEM  code  had  to  be  implemented.  They  define,  in  a  special  format,  coefficients  of  the  volume 
and  boundary  integrals  introduced  in  previous  section.  As  regards  the  contact  condition, 
it  was  assumed  that  contact  can  take  place  at  any  point  of  the  boundary  on  which  static 
boundary  conditions  are  applied.  Therefore,  while  the  integrals  over  this  part  of  boundary 
are  evaluated,  the  program  examines  whether  an  integral  point  is  in  contact  with  the  flat. 
If  penetration  is  greater  than  zero,  then  integrals  corresponding  to  contact  and  friction 
are  added  to  the  coefficients  of  the  stiffness  matrix  and  the  right-hand  side.  Additionally, 
adhesion  integrals  are  being  evaluated  right  outside  the  contact  zone.  In  order  to  properly 
capture  the  strongly  nonlinear  separation-dependent  adhesion  forces,  very  fine  integration 
schemes  are  being  used. 

To  enable  automatic  generation  of  meshes  for  asperity  analysis,  five  additional  programs 
were  prepared.  They  generate  customized  grid  files  for: 

•  3D  axisymmetric  asperity  (cosine  hill), 

•  3D  axisymmetric  asperity  (spherical), 

•  2D  asperity  (cosine  hill), 

•  2D  asperity  (cylindrical), 

•  2D  trapezoidal  asperity. 

The  first  two  programs  make  use  of  an  in-house  GAMMA3D  mesh  generator  .  All  of 
them  provide  generation  of  meshes  with  first  and  second  order  of  geometry  approximation. 
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5  Basic  Verification  of  Numerical  Models 


To  confirm  reliability  of  the  code  and  material  models  several  numerical  tests  were  performed. 
Selected  tests  are  described  in  this  section. 

Several  basic  tests  were  carried  out  to  verify  the  linear  elasticity  formulation.  They  were 
performed  for  both  3D  and  2D  problems  with  second  and  first  order  geometry  approxima¬ 
tion.  The  results  of  tests  were  compared  with  analytical  solutions.  Satisfactory  results  were 
obtained  for  all  the  polynomial  orders  (1  through  8). 

The  objective  of  the  next  group  of  tests  was  to  verify  incremental  formulation  of  the 
viscoplastic  problem.  They  were  carried  out  for  alloy  B1900-fHf  at  temperature  871°C. 
Material  constants  as  well  as  experimental  results  for  this  material  are  given  in  reference 
[10].  For  Bodner-Partom  model  the  material  constants  are  as  follows: 


D0  = 

104  s-1 

m  i  = 

0.270  MPa-1 

n  — 

1.03 

m2  = 

1.52  MPa-1 

Zo  = 

2400  MPa 

ri  = 

r2  =  2 

Z2  = 

2400  MPa 

A\  = 

A2  =  0.0055  s"1 

Zi  = 

3000  MPa 

E  = 

142  GPa 

Z3  = 

1150  MPa 

v  = 

0.0805 

Some  of  these  tests  are  listed  below: 

(a)  Solution  of  the  uniaxial  tension  for  a  purely  elastic  body  using  the  incremental 
formulation.  The  results  were  the  same  as  obtained  with  the  total  formulation. 

(b)  Solution  of  the  uniaxial  creep  test  for  an  elasto-visco-plastic  body.  The  results 
were  compared  with  an  experiment  presented  in  reference  [56],  Certain  discrep¬ 
ancy  of  results  was  observed.  This  discrepancy  was  caused  by  an  erroneous  value 
of  the  Young  modulus  given  in  reference  [56].  After  a  correction  of  Young  modulus 
( E  =  132  GPa )  the  numerical  and  experimental  results  agreed  satisfactorily  (see 
Fig.  5.1).  The  test  verified  mathematical  and  numerical  models  for  this  simple 
loading  (uniaxial  tension). 

(c)  Next  three  tests  were  carried  out  in  order  to  examine  the  performance  of  the 
computer  program  for  more  complicated  loading  histories.  The  tests  included: 

•  Uniaxial  cyclic  tension  and  compression  shown  in  Fig.  5.2. 

•  Uniaxial  loading  for  5,000  s  and  relaxation  for  next  5,000  s.  See  Fig.  5.3. 

•  Loading  for  1,000  s  and  creep  for  next  1,000  s.  See  Fig.  5.4. 

The  above  basic  tests  verified  the  formulation  of  the  viscoplastic  problem  as  well  as  its 
implementation  in  the  code. 
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Figure  5.1:  Experimental  and  numerical  results  of  uniaxial  tension  of  a  specimen  with  cor¬ 
rected  Young  modulus 
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Figure  5.2:  Cyclic  loading  test  of  the  viscoplastic  model  -  10  cycles. 
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Figure  5.3:  Relaxation  test  of  the  viscoplastic  model 
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Figure  5.4:  Creep  test  of  the  viscoplastic  model 
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6  Verification  of  Numerical  Models  of  Asperity 


Numerical  modeling  of  response  of  surface  asperities  to  contact  and  friction  loads  is  one  of 
basic  components  of  the  new  asperity- based  interface  models  developed  in  this  project.  In 
order  to  verify  correctness  of  finite  element  asperity  simulations,  we  performed  several  tests 
and  comparisons,  in  particular: 

•  modeling  of  elastic  sphere  in  contact  with  a  rigid  flat,  (Hertz  problem).  Analytical 
solution  is  available  for  this  problem  [48,91], 

•  numerical  modeling  and  experimental  measurements  for  custom-made  asperities  with 
strongly  pronounced  nonelastic  properties. 

Details  of  these  tests  are  presented  further  in  this  section. 


6.1  Elastic  Sphere  in  Contact  with  a  Rigid  Flat 

In  order  to  verify  the  contact  algorithm  and  the  nonlinear  solution  procedure,  a  finite  element 
solution  was  obtained  for  the  contact  of  elastic  sphere  with  a  rigid  flat.  The  finite  element 
solution  of  this  problem  was  compared  with  theoretical  solution  due  to  Hertz  [48.91].  For 
a  given  sphere  of  radius  R  and  prescribed  normal  displacements  of  the  flat  equal  w ,  the 
theoretical  predictions  of  the  contact  radius  r,  contact  area  A  and  total  load  P  are  given  by: 


r 

A 

P 


VRw 

■kRw 


3(1-^) 


Run? 


The  above  problem  was  solved  numerically  using  the  following  dimensionless  data: 


R  =  1.0 

E  =  1000. 

v  =  0.3 

w  =  .001,  .002,  .005,  .01,  .02 
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Due  to  a  localized  nature  of  the  contact,  the  finite  element  mesh  was  defined  only  for  a 
section  of  the  sphere  around  the  contact  zone.  The  problem  was  solved  using  the  Newton 
procedure  combined  with  adaptive  mesh  refinement.  An  example  of  the  final  refined  mesh 
obtained  for  tu=.02  is  shown  in  figure  6.1  (deformed  configuration  is  displayed  and  only 
boundary  elements  are  shown  for  clarity).  The  same  mesh  with  isosurfaces  of  vertical  dis¬ 
placement  is  presented  in  figure  6.2,  the  slicing  plane,  with  stress  <7yy,  is  shown  in  figure  6.3, 
and  the  error  indicators  projected  on  two  slicing  planes  are  displayed  in  figure  6.4. 

The  results  obtained  numerically  compare  favorably  with  numerical  predictions.  A  de¬ 
tailed  comparison  of  theoretical  and  numerical  results  are  shown  in  table  6.5,  and  the  graph¬ 
ical  comparisons  of  predicted  contact  area  and  total  load  are  shown  in  figure  6.6. 

6.2  Experimental  Studies  of  Models  of  Asperity 

In  order  to  verify  numerical  simulation  of  nonelastic  behavior  of  asperities,  several  experi¬ 
mental  measurements  were  performed  and  then  compared  with  numerical  predictions. 

These  tests  included  simple  tension  and  compression  problems  designed  to  verify  nonelas¬ 
tic  material  constants  for  aluminum,  as  well  as  contact  tests  for  two  types  of  custom-made 
asperities.  In  this  phase  of  experiments,  which  we  refer  to  as  Phase  I,  custom  asperities  were 
chosen  in  order  to  eliminate  random  surface  factor  from  the  comparisons.  In  Phase  II.  which 
will  be  performed  in  the  next  year  of  the  project,  real  random  surfaces  will  be  considered. 

The  objective  of  this  Phase  I  experimental  study  is  to  study  the  deformation  of  contact 
surface  asperities  and  to  verify  the  analytical  prediction  by  experimentation.  In  design  of  the 
experimental  study,  it  is  initially  conceived  that  the  test  is  to  be  carried  out  under  a  small 
normal  load  (500  Lbf)  condition.  A  controlled  surface  asperities  are  machined  onto  both 
surfaces  of  an  aluminum  block.  In  the  test  arrangement,  the  aluminum  block  is  sandwiched 
between  two  hardened  steel  blocks  with  smooth  surfaces,  and  the  steel-aluminum-steel  block 
assembly  is  compressibly  loaded  to  approximately  500  Lbf.  The  deformation  of  asperities 
on  the  aluminum  block  is  monitored  during  loading.  After  completing  the  normal  loading, 
a  horizontal  load  is  slowly  applied  to  the  aluminum  block  until  the  block  begins  to  slip. 
An  estimation  of  the  coefficient  of  friction  can  thus  be  made  by  using  the  measured  normal 
and  horizontal  loads,  and  be  compared  with  that  from  the  analytical  prediction.  A  test 
apparatus  is  built  for  this  purpose,  and  tests  are  carried  out  with  this  apparatus.  After 
reviewing  the  test  results,  it  is  decided  that  the  deformation  of  surface  asperities  on  the 
aluminum  specimen  is  too  large  for  the  purpose  of  model  verification.  The  shape  of  asperity 
is  changed,  and  the  tests  are  carried  out  under  a  normal  load  of  approximately  2,000  Lbf. 
The  results  from  both  arrangements  are  reported  in  this  section. 
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ELASTIC  SPHERE  IN  CONTACT  WITH  A  RIGID  FLAT 


E  =  1.0000E+03  nu=  0.3000  R-  1.0000 


THEORY 


Displ  Cont .  rad.  Area  Load  Max.  pres. 


0.00200 

0.00500 

0.01000 

0.02000 

0.04472 

0.07071 

0.10000 

0 . 14142 

0.006283 

0.015708 

0.031416 

0.062832 

1 . 3105E-01 

5 . 1803E-01 
1.4652E+00 

4 . 1442E+00 

3.1286E+01 
4 . 9468E+01 

6 . 9958E+01 
9.8936E+01 

NUMERICAL 

Displ 

Cont.  rad. 

Area 

Load 

Max.  pres. 

0.00200 

0.00500 

0.01000 

0.02000 

0.04237 

0.07001 

0.10092 

0.14809 

0.00564 

0.0154 

0.0320 

0.0689 

1 . 5400E-01 

5 . 5400E-01 

1 . 5400E+00 

4 . 5000E+00 

4 . 3600E+01 
5 . 4000E+01 
7 . 9200E+01 
1 . 0840E+02 

Figure  6.5:  Comparison  of  theoretical  and  numerical  results  for  the  Hertz  problem 
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a) 


b) 


Figure  6.6:  Comparison  of  theoretical  and  numerical  results  for  the  Hertz  problem;  (a) 
displacement  versus  contact  area,  and  (b)  displacement  versus  contact  load 
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6.2.1  The  Test  Apparatus 


A  sketch  of  the  apparatus  is  shown  in  Fig.  6.7.  An  aluminum  (6061,  T4)  block  of  1’  x  1’  x 
0.25”  is  sandwiched  between  two  steel  blocks  of  the  same  dimension  as  shown  in  the  figure. 
The  V-shaped  grooves  of  angle  45  deg,  pitch  spacing  of  0.1  inch,  and  depth  of  0.125  inch 
are  machined  with  a  specially  designed  cutter  on  both  surfaces  of  the  aluminum  block  to 
simulate  the  surface  asperities.  A  sketch  of  the  grooves  is  shown  in  Fig.  6.8.  The  surfaces  of 
the  steel  blocks  are  machined  smooth  and  (water)  quench  hardened  to  RC-30.  It  should  be 
mentioned  that,  due  to  the  angle  of  the  cutting  tool  and  the  pitch  spacing,  the  tips  of  the 
grooves  are  not  in  a  plane,  the  heights  of  tips  vary  alternatively  as  shown  in  the  photograph 
(Fig.  6.8)  to  be  discussed  in  a  later  section. 

A  normal  load  is  applied  to  the  specimen  assembly  through  a  mechanical  screw  jack  from 
the  bottom  of  the  apparatus.  The  normal  load  is  monitored  with  a  small  “load  transducer” 
of  capacity  500  Lbf.  The  normal  deformation  of  the  simulated  asperities  on  aluminum  block 
is  measured  with  a  “proximate  sensor”  with  an  operation  range  from  0  to  0.1  inch.  The  sand¬ 
wiched  specimen  assembly  is  first  installed  in  position,  and  an  initial  load  of  approximately 
50  Lbfs  is  applied  to  the  specimen  prior  to  “zero  adjustment”  of  the  recording  instrument 
(an  X-Y  plotter).  A  total  normal  load  of  500  Lbs  is  then  applied  to  the  specimen  at  a  rate  of 
approximately  10  Lbf  per  minute.  After  the  normal  load  has  reached  500  Lbf,  a  horizontal 
load  (by  lead  blocks  and  beads)  is  slowly  applied  to  the  aluminum  block  until  the  block  starts 
to  slide.  The  horizontal  load  is  monitored  with  a  ring-shaped  load  transducer  (laboratory 
built)  as  shown  in  Fig.  6.7. 

Note  that  even  with  very  precise  calibration  of  the  test  apparatus,  certain  compliance 
or  “setting  in”  occurs  during  loading  and  pollutes  the  measurements.  This  is  especially 
true  in  case  of  high  loads  and  very  small  displacements  considered  in  this  experiment.  It  is 
a  standard  practice  in  experimental  tests  to  discard  the  initial  part  of  the  load  curve  and 
appropriately  translate  the  remaining  part.  In  this  section,  we  present  results  in  the  “raw” 
form,  but  numerical  comparisons  refer  to  corrected  graphs. 

6.2.2  Tests  Results  From  the  Above  Apparatus 

Two  tests  are  carried  out.  A  representative  force- deformation  plot  is  shown  in  Fig.  6.9.  It 
is  seen  that  there  is  an  appreciable  amount  of  plastic  deformation  of  the  surface  asperities 
on  the  aluminum  block  (neglect  the  bulk  deformation  of  both  aluminum  and  steel  blocks) 
when  the  block  assembly  is  loaded  to  500  Lbf.  A  horizontal  force  is  then  slowly  applied  to 
the  aluminum  block.  It  is  visually  observed  that  the  aluminum  block  starts  to  slip  when  the 
horizontal  load  reaches  108  Lbf.  Since  there  are  two  contact  surfaces  between  the  aluminum 
block  (with  asperities  on  both  surfaces  as  shown  in  Fig.  6.8)  and  steel  blocks  (with  smooth 
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surfaces)  in  the  test  arrangements,  the  nominal  coefficient  of  friction  is  (108  /2)/500  =  0.108. 

6.2.3  Deformation  of  Asperities  Under  a  Larger  Normal  Load 

The  asperities  shown  in  Fig.  6.8  have  pointed  tips,  the  deformation  of  tip  is  difficult  to 
calculate.  It  is  then  decided  to  change  the  geometry  of  asperity  as  those  shown  in  Figs. 
6.8(b)  and  6.8  (c).  In  Fig.  6.8  (b),  the  asperity  is  modeled  as  a  45-deg  grooved  with  a 
truncated  tip.  In  Fig.  6.8  (c),  the  asperity  is  modeled  by  spaced  circular  rods.  The  force- 
deformation  relationship  for  the  specimen  with  asperities  as  shown  in  Fig.  6.8  (b)  and  6.8 
(c)  are  measured.  The  test  arrangements  are  the  same  as  that  described  in  the  previous 
section.  Since  the  maximum  load  in  these  test  are  much  higher  than  500  Lbf,  the  tests  are 
carried  out  and  the  representative  force-deformation  curves  are  shown  in  Figs.  6.10  and 

6.11,  respectively.  A  photograph  of  the  cross-section  of  the  deformed  grooves  is  shown  if  Fig. 

6.12.  It  is  seen  that  only  the  alternate  grooves  were  deformed.  There  are  seven  grooves  on 
each  surface  of  the  aluminum  block,  therefore  only  eight  (four  on  each  surface)  of  them  are 
deformed.  Since  the  maximum  horizontal  load  required  to  move  the  aluminum  block  in  this 
exceeds  the  capacity  of  the  ring  load  cell  used  in  previous  section,  the  horizontal  pulling  test 
has  yet  to  be  carried  out  (we  need  to  build  a  new  load  cell  and  loading  frame). 

Finally,  a  test  for  the  tensile  property  of  the  aluminum  used  in  this  study  is  carried  out. 
The  stress-strain  curve  from  an  aluminum  (6061,  T4)  specimen  is  shown  in  Fig.  6.13.  Again, 
the  test  is  carried  out  at  a  slow  loading  rate  (approximately  20  Lbf  per  minute). 


88 


Figure  6.7:  A  sketch  of  test  apparatus  (not  to  scale). 
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Figure  6.8:  Different  models  of  asperities  studied  experimentally. 
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Figure  6.9:  Load-displacement  cun 
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Figure  6.12:  A  photograph  of  deformed  asperity 

6.3  Numerical  Simulation  of  Experimental  Measurements 

6.3.1  Viscoplastic  Uniaxial  Stress  State 

The  aim  of  this  final  test  was  to  determine  viscoplastic  material  constants  for  aluminum  6061 
T4.  Bodner-Partom  model  of  viscoplastic  materials  uses  a  total  of  14  material  constants. 
However,  for  aluminum  alloys  at  room  temperature  only  7  material  constants  are  of  primary 
importance.  Their  values  obtained  from  reference  [10]  for  a  different  heat  treatment  (T6) 
are  listed  below: 

E  =  73.9  GPa 

v  =  0.33 

zq  =  450  MPa 

z\  —  550  MPa 

mi  =  0.12  MPa-1 

D  =  108s_1 

n  =  5.0 

In  order  to  verify  these  values,  a  tension  test  was  performed  and  compared  with  numerical 
results.  The  test  indicated  that  the  6061  T4  sample  has  slightly  different  values  of  Young 

modulus  and  the  kinematic  hardening  parameter.  These  values  are: 

E  =  65.0  GPa 

n  =  5.8 

These  material  constants  were  used  in  further  computations.  Fig.  6.14  shows  the  plot 
of  stress-strain  relation  for  initial  and  modified  material  constants  in  comparison  with  the 
experimental  results. 
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6.3.2  Viscoplastic  Cylindrical  Asperity 

Aluminum  cylinder  was  used  as  a  model  of  a  ”  macro  asperity”.  Numerical  and  experimental 
tests  were  carried  out.  It  was  assumed  that  the  cylinder  can  be  analyzed  as  a  2D  case.  The 
original  mesh  used  for  discretization  of  a  section  of  a  circle  is  shown  in  Fig.  6.15.  After 
the  first  solution  pass  the  mesh  whs  automatically  refined.  It  is  shown  in  Fig.  6.16.  The 
viscoplastic  results  obtained  at  both  meshes  as  well  as  a  purely  elastic  solution  for  the  first 
mesh  are  compared  in  Fig.  6.17.  The  conclusions  are  that: 

•  behavior  of  the  specimen  under  applied  loading  is  almost  elastic, 

•  viscoplastic  solution  is  reasonable  -  it  gives  smaller  values  of  the  contact  force, 

•  the  refined  mesh  gives  results  which  can  be  treated  as  a  final  numerical  solution  (the 
difference  between  solutions  obtained  at  coarse  and  fine  meshes  is  small.) 

The  numerical  and  experimental  solutions  are  compared  in  Fig.  6.18.  Note  that  the 
experimental  results  have  been  rescaled.  Instead  of  the  total  force  -  total  displacement 
relation  measured  in  the  experiment,  Fig.  6.18  shows  force  per  unit  length  of  upper  half- 
cylinder.  Moreover,  metric  units  were  used. 

It  can  be  observed  that  the  numerical  results  compare  favorably  with  experimental  mea¬ 
surements.  Recall,  however,  that  the  experimental  data  were  translated,  to  correct  for 
settling  in  the  apparatus  (see  Section  6.1). 

6.3.3  Viscoplastic  Custom  Surface  Model 

Comparisons  of  results  for  a  model  truncated  V-shaped  of  asperity  are  described  in  this 
section.  Two  meshes  which  were  used  for  discretization  of  this  specimen  are  shown  in  Fig. 
6.19  and  6.20.  Numerical  results  for  both  meshes  are  shown  in  Fig.  6.22  and  6.21.  While  for 
the  cylindrical  model  the  behavior  of  the  material  was  almost  elastic,  in  this  case  the  influence 
of  yielding  was  significant.  Maximum  nonlinear  strain  was  about  80%.  It  means  that  at 
least  locally  solution  is  beyond  the  theory  of  small  strains  which  we  use.  The  comparison  of 
numerical  results  and  experimental  ones  is  shown  in  Fig.  6.23.  A  good  agreement  of  results 
can  be  observed  for  both  models  of  asperity. 

Some  of  unevitable  sources  of  discrepancy  are  the  following: 

•  error  of  experimented  measurements, 

•  modeling  of  inelastic  behavior  of  the  material  beyond  the  small  deformation  range, 
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•  unknown  deformation  history  of  the  specimens,  prior  to  the  experiment, 

•  errors  of  time-integration  and  other  numerical  errors. 

Comparison  of  numerical  and  experimental  results  as  well  as  the  basic  numerical  tests 
lead  to  a  conclusion  that  the  adaptive  finite  elements  represents  nonelastic  deformation  of 
asperities  with  sufficient  accuracy. 
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Figure  6.14:  Comparison  of  stress-strain  relation  behavior  for  uniaxial  tension  test. 
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Figure  6.15:  Original  discretization  of  the  cylindrical  sample 


Figure  6.16:  Refined  discretization  of  the  cylindrical  sample 


Figure  6.17:  Numerical  results  for  the  cylinder 


Figure  6.18:  Comparison  of  numerical  and  experimental  results  for  the  cylindrical  sample 
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|  Figure  6.20:  Refined  discretization  of  the  second  specimen 
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Figure  6.22:  Numerical  results  for  the  second  specimen  obtained  on  two  different  meshes 
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7  Studies  of  Asperity- Based  Models  of  Contact  and 
Friction 

In  this  section  we  present  examples  of  constitutive  interface  models,  developed  through  a 
complete  statistical  homogenization  procedure.  The  basic  questions  addressed  here,  besides 
the  quality  of  these  models,  are  as  follows: 

•  what  is  the  dominant  type  of  asperity  deformation  (elastic  or  nonelastic), 

•  what  is  the  effect  of  adhesion  at  the  asperity  level  and  after  homogenization, 

•  what  is  the  effect  of  surface  roughness, 

•  how  do  theoretical  predictions  compare  with  experimental  measurements. 

To  answer  these  questions  four  different  types  of  surfaces  were  analyzed.  The  first  three 
were  made  of  aluminum.  Their  profiles  were  taken  from  literature  [44].  The  fourth  type  of 
surface  was  made  of  steel.  The  pressure-approach  relationship  was  measured  experimentally 
and  compared  with  its  numerical  prediction.  The  profile  was  scanned  with  an  electronic 
microscope.  Having  the  profile  of  a  surface  as  a  discrete  function  it  was  possible  to  compute 
its  statistical  characteristics  such  as  medium  height,  standard  deviations  of  height,  slope 
and  curvature.  Those  statistical  characteristics  were  used  to  compute  expected  values  of 
real  area  and  force  of  contact  as  functions  of  approach.  Deterministic  solutions  to  a  contact 
of  two  asperities  were  computed  by  the  hp  adaptive  code  using  different  assumptions  about 
shape  of  asperities,  its  material  model  and  with  or  without  taking  into  account  the  adhesion 
forces.  The  results  include  also  estimation  of  friction  coefficient. 

7.1  Simulation  of  a  Greenwood- Williamson  Asperity-Based 
Contact  Model 

The  results  of  elastic  solution  of  spherical  contact  problem  were  used  for  initial  tests  of 
a  complete  statistical  homogenization  procedure.  In  particular,  we  compared  our  finite 
element  based  predictions  with  a  classical  asperity-based  contact  model  due  to  Greenwood 
and  Williamson  [44].  This  is  one  of  historically  first  asperity-based  models,  derived  under 
the  following  simplifying  assumptions: 

•  the  tips  of  asperities  are  spherical, 

•  all  asperities  have  the  same  radius  of  curvature  R, 
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•  asperity  height  is  a  random  variable  with  Gaussian  distribution, 

•  contact  is  elastic  and  described  by  the  theoretical  solution  of  the  Hertz  problem. 

Although  this  model  is  much  simpler  than  the  ultimate  objective  of  this  project  (random 
asperity  height  and  curvature,  nonelastic  deformation,  etc.),  it  provides  a  very  good  initial 
test  for  the  homogenization  procedure.  Importantly,  our  modeling  of  this  problem  differs 
from  Greenwood’s  approach  in  that: 

•  solution  of  elastic  contact  problem  has  been  obtained  from  finite  element  modeling, 
rather  than  from  Hertz  theory, 

•  expected  values  of  contact  load  and  area  have  been  calculated  using  general  numerical 
quadrature,  instead  of  analytical  integration. 

The  simulation  of  the  Greenwood- Williamson  model  was  performed  for  the  following  set 
of  parameters: 

R  —  0.1414  mm 
av  =  7.07106  x  10-4  mm 
v  =  0.3 

E  =  321.7335  kG\mm 2 
n  =  300  peaks\mm 2 

where  R  is  the  radius  of  asperity  tips,  crp  is  the  standard  deviation  of  asperity  height  and  n 
is  the  surface  density  of  asperity  peaks.  The  results  are  presented  in  figure  7.1,  which  shows 
respective  comparisons  of  load-separation  and  load  vs. contact  area  curves.  It  can  be  seen 
that  the  difference  introduced  by  numerical  modeling  of  asperity  and  numerical  integration  of 
expectation  values  is  within  acceptable  bounds  (note  that  these  differences  could  be  further 
reduced  by  application  of  finer  meshes  for  the  asperity  contact  problem.) 

7.2  Effects  of  Asperity  Shape 

Next  test  of  the  homogenization  procedure  was  performed  for  an  aluminum  interface  with 
surface  characterized  by  the  following  parameters:  a  —  1.3  pm,  &  =  0.13,  a  =  0.018  pm~l, 
Dp  =  0.00044  peaks  /pm2.  During  this  test  the  force  and  area  of  contact  for  a  single  asperity 
were  computed  in  two  ways,  using 

•  homogenization  procedure  based  on  elastic  solution  for  sperical  asperities,  with  random 
distribution  not  only  of  heights  but  also  of  curvatures,  (Hertz  solution  was  used). 
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Figure  7.1:  Comparison  of  Greenwood- Williamson  theory  with  numerical  predictions;  (a) 
load- separation  curve  and  (b)  load- area  curve. 


•  Homogenization  procedure  based  on  numerical  solutions  for  elastic,  sinusoidal  asperi¬ 
ties.  These  asperities  better  simulate  the  actual  shape  of  a  rough  surface. 

Numerical  solution  for  sinusoidal  asperities  was  computed  with  taking  advantage  of  the 
fact  that  the  problem  is  axisymmetric.  Comparison  of  the  computed  contact  force  and 
axea  versus  penetration  is  shown  in  figures  7.2  and  7.3.  Results  for  contact  of  two  surfaces 
obtained  by  those  two  models  are  compared  in  figure  7.4. 

The  results  of  this  test  show  that  influence  of  asperity  geometry  (sinusoidal  rather  of 
spherical)  is  important,  especially  for  higher  loadings.  However,  after  homogenization,  these 
differences  axe  less  pronounced.  This  is  because  major  contribution  to  expected  values  of 
contact  pressure  and  area  come  from  asperities  deformed  at  their  tips  only. 

7.3  A  very  smooth  engineering  surface 

The  next  case  analyzed  is  normal  contact  of  two  very  smooth  aluminum  surfaces,  with  profiles 
corresponding  to  that  shown  by  Greenwood  and  Williamson  in  Reference  [18],  Figure  5.  The 
material  is  aluminum  alloy  (6061,  T4),  with  constitutive  parameters  defined  in  Section  6. 
The  surface  roughness  corresponds  approximately  to  a  very  well  polished  bearing  surface, 
and  is  defined  by  the  following  parameters: 

a  —  0.013  fj.m 
<7  =  0.013 
a  =  0.018  fim~l 
Dp  =  0.058  p.m~2 

To  represent  relative  surface  roughness,  it  is  convenient  to  use  a  so-called  plasticity  index, 
defined  by  Greenwood  and  Williamson  [44]  as: 

where  E'  =  E(  1  —  i/2),  H  is  the  material  hardness,  and  R  is  the  asperity  radius.  According 
to  arguments  presented  by  Greenwood  and  Williamson,  for  the  values  of  *P  smaller  than 
1  the  surface  deformation  is  essentially  elastic  at  a  wide  range  of  normal  loads,  while  for 
$  >  1  plastic  deformation  can  be  expected.  Since  in  the  present  approach  the  asperity 
peak  curvature  is  a  random  variable,  we  obtain  a  range  of  plasticity  indices  corresponding  to 
various  asperities  under  consideration.  For  the  surface  considered  here,  the  value  of  plasticity 
index  varies  between  'P  =  0.8  and  »P  =  3.5. 
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Figure  7.3:  Area  of  contact  for  a  sinusoidal  asperity 
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Figure  7.4:  Nominal  pressure  vs.  approach  for  an  interface;  spherical  (a)  and  sinusoidal  (b) 
asperities 


In  order  to  study  the  effect  of  adhesion  for  smooth  surfaces,  considered  cases  with  and 
without  adhesion  forces  (clean  and  contaminated  surfaces,  respectively).  The  model  of  ad¬ 
hesion  corresponds  to  the  one  discussed  in  this  work,  with  the  surface  energy  of  adhesion 
A7  =  |J/m2[52]. 

The  first  step  of  the  analysis  was  to  use  finite  element  method  to  model  loading  of 
severed  individual  asperities,  with  radii  spanning  the  effective  support  of  the  probability 
density  function,  namely  from  R=15  microns  to  R=300  microns.  A  typical  mesh  for  such 
an  analysis  is  shown  in  Figure  7.5.  A  very  fine  refinement  around  the  perimeter  of  the 
contact  zone  was  needed  to  correctly  resolve  the  strongly  nonlinear  adhesion  forces.  A 
sample  solution  for  a  selected  asperity  of  radius  R=15.4  microns  is  shown  in  the  form  of 
load-deflection  curves  in  figure  7.6.  It  is  of  interest  to  notice,  that  in  presence  of  adhesion 
(clean  surfaces  in  vacuum),  certain  attractive  force  develops  before  the  asperities  come  into 
contact.  In  these  computations,  elastic  asperity  bahavior  was  assumed. 

After  sending  the  finite  element  results  through  the  statistical  homogenization  package, 
one  obtains  pressure- approach  and  contact  area  approach  curves  as  shown  in  Figures  7.7  and 
7.8.  Interestingly,  due  to  a  specific  statistical  distribution  of  asperity  heights,  the  effect  of 
adhesion  is  more  pronounced  on  the  surface  level  than  for  a  single  asperity  and,  for  clean 
smooth  surfaces  in  vacuum,  a  considerable  attractive  force  can  be  expected. 
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Figure  7.6:  Smooth  aluminum  surface:  load  deflection  curves  for  a  single  asperity 
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Figure  7.7:  Smooth  aluminum  surface:  normal  stress  versus  approach  curve 

One  comment  is  in  place  regarding  the  “approach”  axis  in  Figure  7.6.  By  definition, 
the  approach  between  two  surfaces  is  measured  from  the  moment  when  the  surfaces  (the 
highest  asperity  tips)  come  into  contact.  For  real  surfaces  this  initial  contact  is  relatively 
easy  to  capture.  However,  in  theoretical  analysis,  with  Gaussian  distribution  of  surface 
height,  it  is  not  the  case  -  the  maximum  asperity  height  is  not  clearly  defined  (in  fact  it 
grows  to  infinity,  with  probability  density  vanishing  to  zero).  As  a  consequence,  the  selection 
of  the  zero  point  on  the  “approach”  axis  is  somewhat  arbitrary  -  for  practical  purposes  it 
may  be  chosen  as  a  point  where  some  noticeable  normal  load  develops  on  the  interface. 
In  Figure  7.7,  for  illustration  purposes,  the  approach  was  measured  from  a  rather  large 
separation  between  the  surfaces.  Importantly,  this  approach  measurement  does  not  affect 
the  relationship  between  the  real  contact  area,  normal  load  and  normal  stiffness,  which  is 
of  greater  practical  importance.  For  the  surface  under  consideration  the  real  contact  area 
calculated  for  increasing  load  is  shown  in  Figure  7.9. 

As  mentioned,  the  above  asperity  analysis  was  performed  using  small  deformation  theory, 
with  additional  assumption  of  elastic  asperity  deformation.  Therefore  a  question  arises:  Of 
all  the  asperities  in  contact,  how  many  satisfy  the  above  assumptions?  An  answer  to  this 
question  is  shown  in  Figure  7.10,  which  presents,  for  increasing  approach,  the  percentage  of 
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Figure  7.9:  Smooth  aluminum  surface:  contact  area  for  increasing  normal  pressure 

asperities  that  are  in  purely  elastic  contact,  as  well  as  the  percentage  of  asperities  within 
the  range  of  small  deformation  theory  (|e|max  <  8%).  Clearly,  the  surface  deformation 
remains  well  within  the  range  of  small  deformation  theory.  However,  the  number  of  nonelastic 
asperities  increases  with  increasing  approach  (load).  This  corresponds  quite  well  to  the 
behavior  suggested  by  the  range  of  plasticity  index  for  the  surface  under  consideration.  In 
order  to  assess  the  importance  of  nonelastic  effects,  an  additional,  fully  nonlinear  derivation 
was  performed,  for  the  case  without  adhesion.  A  comparison  of  approach-pressure  curves 
obtained  with  purely  elastic  and  elasto-plastic  approaches  is  shown  in  figure  7.10.  The 
maximum  difference  between  elastic  and  elasto-plastic  results  was  less  that  5%. 

It  can  be  noted,  that  a  very  localized  occurrence  of  plastic  deformation  in  an  asperity 
does  not  strongly  affect  the  values  of  load  and  contact  area,  and  that  the  load-deflection  and 
load-area  curves  obtained  with  the  elastic  theory  are  quite  reliable,  even  with  50  percent  of 
asperities  outside  the  purely  elastic  range.  This  is  because  only  for  about  7%  of  asperities 
the  elasto-plastic  solution  was  significantly  different  than  the  purely  elastic  one. 
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itage  of  asperities  satisfying  assumptions  of 


7.4  Studies  of  a  rough  surface 


A  second  example  considered  here  is  that  of  a  hypothetical  rough  engineering  surface.  Sim¬ 
ilarly  as  in  the  previous  case,  the  material  is  aluminum  alloy,  and  the  surface  statistics  is 
defined  by: 


a  —  1.3  pm 
<7  =  0.13 
a  =  0.018  pm-1 
Dp  =  0.00044  pm-2 

For  this  surface,  the  plasticity  index  is  well  over  the  value  of  30.0,  so  that  extensive  plastic 
deformation  can  be  expected.  In  order  to  derive  asperity- based  model  of  the  interface,  several 
asperities  with  various  peak  curvatures  were  analyzed  by  the  finite  element  method.  The 
deformation  was  fully  elasto-viscoplastic,  and  the  asperities  were  subjected  to  normal  load 
executed  by  a  rigid  flat  surface.  Importantly,  for  the  rough  surface  the  results  were  not 
strongly  affected  by  adhesion.  This  effect  corresponds  well  to  the  observations  of  Chang,  et. 
al.  [26-28],  that  the  effect  of  adhesion,  even  for  clean  surfaces,  diminishes  with  increasing 
surface  roughness.  On  the  other  hand,  for  rough  surfaces  the  effect  of  nonelastic  deformation 
is  very  significant  -  this  can  be  seen  from  the  comparison  of  elastic  and  elasto-plastic  pressure- 
approach  curves  shown  in  Figure  7.11.  Clearly,  the  difference  between  elastic  and  elasto- 
plastic  solution  is  very  pronounced.  A  similar  effect  can  be  observed  on  the  plot  showing 
the  dependence  of  real  contact  area  on  the  normal  pressure  -  Figure  7.13.  Importantly,  this 
difference  will  have  a  very  strong  effect  on  the  values  of  the  coefficient  of  friction. 

It  is  of  interest  to  look  again  at  the  percentage  of  asperities  within  elastic  range  and  within 
small  deformation  range.  These  numbers  are  shown  in  figure  7.14.  It  can  be  observed,  that  for 
rough  surfaces  even  a  very  small  load  causes  significant  plastic  deformation  -  the  percentage 
of  elastic  asperities  is  essentially  zero.  Moreover,  at  higher  loads,  some  asperities  experience 
strain  levels  exceeding  the  range  of  small  strain  theory  -  an  important  observation,  especially 
for  more  flexible  polymer  surfaces. 
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Figure  7.13:  Rough  aluminum  surface:  real  contact  area  for  increasing  normal  pressure 
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8  Experimental  Verification  of  Asperity-Based  Con¬ 
tact  Models 

The  ultimate  test  of  new  asperity-based  models  of  contact  interfaces  comes  in  applications 
to  predicting  phenomenological  behavior  of  real  engineering  surfaces. 

To  verify  these  new  models,  a  special  experiment  was  designed  and  performed.  The  main 
objective  was  to  compare  experimental  force- deflection  relationship  with  the  asperity- based 
theoretical  prediction.  Additionally,  measured  and  computed  coefficients  of  friction  were  also 
compared.  The  following  subsections  describe  the  experiment  and  numerical  computations. 

8.1  Experimental  Samples,  Apparatus,  and  Measurements 

A  major  difficulty  in  experimental  verifications  of  models  of  frictional  interfaces  is  caused  by 
relatively  small  height  of  surface  asperities  and  resulting  small  compliance  of  the  interface. 
Under  these  conditions,  it  is  rather  difficult  to  avoid  pollution  of  the  results  by  such  effects 
as: 

•  deformation  of  a  bulk  material, 

•  departure  of  experimental  contacting  surfaces  from  ideally  flat,  which  spoils  homoge¬ 
neous  distribution  of  contacting  asperities  and  creates  additional,  plate-bending  type 
compliance, 

•  compliance  of  the  loading  apparatus, 

•  departure  of  loads  from  purely  axial,  etc. 

In  order  to  minimize  these  negative  effects,  special  specimens  and  apparatus  were  used, 
as  described  below. 

8.1.1  Specimen  Preparation  and  Experimental  Arrangement: 

The  disk  specimens  of  1/8  in.  thickness  and  1  in.  diameter  were  cut  from  a  cold  rolled  1020 
steel  rod.  The  rod  has  a  hardness  of  89  /?&,  and  the  tensile  stress-strain  relationship  of  the 
rod  material  is  shown  in  figure  8.1.  The  material  has  an  elastic  modulus  E  =  30  x  106  psi, 
and  a  yield  stress  approximately  42  x  103  psi. 

Ten  disks  were  used  in  each  test.  The  disks  were  first  lapped  with  a  milling  machine  to 
ensure  the  same  surface  flatness  among  them.  Both  surfaces  of  the  disks  were  then  sand¬ 
blasted  with  steel  and  glass  beads  to  produce  an  artificial  surface  roughness.  In  the  test 
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arrangement,  ten  disks  were  stacked  to  form  a  column  of  nine  (9)  contact  surfaces  as  shown 
in  figure  8.2.  Note  that,  for  the  first  disk  (no.  1)  and  the  last  disks  (no.  10),  only  the  contact 
surfaces  with  the  adjacent  disks  were  sand-blasted,  the  contact  surfaces  of  these  disks  (no. 
1  and  10)  with  anvils  of  the  testing  machine  were  kept  smooth  as  lapped.  The  purpose  of 
using  the  stacked  column  was  to  increase  the  deformation  of  contact  surfaces  nine  times  so 
that  the  average  value  of  deformation  of  the  contact  surface  asperities  could  be  measured 
with  a  reasonable  accuracy.  As  portrayed  in  the  figure,  the  applied  load  was  measured  with 
a  load  cell,  and  the  displacement  was  measured  with  a  LVDT  which  has  a  resolution  of  0.001 
inch. 


8.1.2  Measurement  of  Surface  Roughness 

The  average  depth  of  surface  asperities  produced  by  sand-blasting  was  first  estimated  with 
a  comparator  which  is  a  device  commonly  used  by  the  sand-blasting  industries.  It  was 
estimated  that  the  surface  which  was  sandblasted  by  the  glass  beads  has  an  average  asperity 
depth  of  20  yu,  and  the  surface  which  was  sandblasted  by  the  steel  beads  has  an  average 
asperity  depth  of  40  /. i .  These  asperities  are  too  small  to  be  monitored  by  a  standard  stylus 
profilometer.  That  is  why  the  surface  roughness  of  the  specimen  were  determined  with  a 
scanning  electronic  microscope  (courtesy  of  NASA-JSC).  The  surface  of  the  disk  specimen 
was  first  scanned  to  determine  the  distribution  of  asperities  (i.e.,  a  top  view).  The  disk 
specimen  was  then  sectioned,  lapped,  polished,  and  scanned  from  the  side  to  determine  the 
depth  variation  of  surface  asperities  (i.e.  a  side  view  profile).  The  results  are  presented  in 
the  following  paragraphs: 

1.  The  top  and  side  view  pictures  of  surface  asperities  produced  by  sand-blasting  the 
surface  with  glass  beads  are  shown  in  Figs.  8.3  and  8.4,  respectively.  The  variation  of 
asperity  depth,  traced  from  Fig.  8.4,  is  shown  in  Fig.  8.5.  In  view  of  Fig.  8.3,  it  is  seen 
that  the  sand-blasting  had  indeed  produced  a  reasonably  homogeneous  distribution  of 
surface  asperities.  From  Figs.  8.4  and  8.5,  it  is  seen  that  the  surface  damage  produced 
by  glass  beads  blasting  is  not  severe;  the  depth  variation  of  asperities  is  gentle  and  their 
distribution  appears  to  be  random.  Based  on  Fig.  8.5,  the  average  depth  of  asperities 
is  estimated  to  be  approximately  15  (i  which  is  slightly  smaller  than  that  estimated  by 
using  the  comparator. 

2.  The  corresponding  information  on  surface  asperities  produced  by  sand-blasting  the 
surface  with  steel  beads  are  shown  in  Fig.  8.6,  8.7,  and  8.8  respectively.  It  is  seen  that, 
although  the  distribution  of  asperities  appears  to  be  uniform,  blasting  the  surface 
with  steel  beads  has  produced  severe  micro-scale  damage  on  the  surface.  The  surface 
asperities  are  more  rough  and  abrupt  than  that  produced  by  sandblasting  the  surface 
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with  glass  beads.  Based  on  Fig.  8.8,  the  average  depth  of  asperities  is  estimated  to  be 
40  n  -  the  same  as  that  from  the  comparator.  Figure  8.9  shows  a  steel  bead  wedged 
between  asperities. 


8.1.3  Experimental  Measurements  of  Contact  Compliance 

The  sand-blasted  disks  were  stacked  together  to  form  a  cylinder,  and  the  cylindrical  specimen 
was  then  installed  to  a  MTS  testing  machine.  A  compressive  force  was  applied  to  the  stacked 
cylinder  at  an  approximate  rate  of  100  Lbf/min.  The  overall  deformation  (between  the 
smooth  surfaces  of  no.  1  and  no.  10  disks)  was  measured  with  a  LVDT  displacement  sensor 
(with  a  resolution  of  0.001  in)  as  shown  in  Fig.  8.2.  Two  tests  for  each  stacked  cylinders, 
which  were  made  of  disks  with  a  glass-bead  blasted  surface  and  with  a  steel-bead  blasted 
surface,  were  carried  out  to  ensure  the  repeatability  of  the  testing.  The  corresponding  results 
are  presented  in  Figs.  8.10  and  8.11,  respectively.  The  following  observation  is  made: 


1.  Figure  8.10  shows  the  compressive  force-displacement  curve  of  a  cylinder  made  by 
stacking  ten  disks  (b/\L  1020)  whose  surfaces  were  sandblasted  with  glass  beads.  The 
cylinder  was  corup  -ssively  loaded  to  1,000  Lbf,  unloaded,  and  then  reloaded  to  3.000 
Lbf.  Note  th'^t  the  diameter  of  the  cylinder  is  1  inch,  the  nominal  compressive  stress 
in  the  cylinder  under  a  load  of  3,000  Lbf  is  approximately  3820  psi  -  well  within  the 
elastic  limit  of  the  cylinder  material.  Furthermore,  each  disk  has  a  thickness  of  0.125 
in.  the  elastic  deformation  of  the  equivalent  solid  cylinder  is  expected  to  be 
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3000  x  (10  x  0.125) 
n  x(0.5)2  x  30  x  106 


=  0.000159  in. 


Referring  to  the  force  deformation  curve  shown  in  Fig.  8.10,  it  is  interesting  to  note 
that  the  initial  portion  of  the  curve  is  nonlinear,  and  there  is  a  permanent  deformation 
of  approximately  0.0007  inch  after  unloading,  and  the  reloading  curve  joins  the  initial 
loading  curve  as  shown  in  the  figure.  At  a  load  of  3,000  Lbf,  the  measured  deformation 
of  the  stacked  cylinder  is  approximately  0.004  inch  which  is  substantially  larger  than  the 
above  calculated  elastic  deformation.  This  phenomenon  implies  that  the  deformation  of 
asperities  between  the  contact  surfaces  (there  are  nine  contact  surfaces  in  the  stacked 
cylinder)  in  the  stacked  cylinder  are  plastic,  and  the  contact  surface  compliance  is 
nonlinear. 


2.  The  corresponding  force-displacement  curve  for  a  cylinder  made  by  stacking  ten  disks 
which  surfaces  are  sandblasted  with  steel  beads  is  shown  in  Fig.  8.11.  It  is  seen  that 
the  characteristics  of  the  force-deformation  relation  is  similar  to  that  displayed  in  Fig. 
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Figure  8.1:  Tensile  stress-strain  relationship  of  a  SAE  1020  steel 

8.10.  The  permanent  deformation  of  the  cylinder  when  unloaded  from  1,000  Lbf  is 
approximately  0.0014  in  -  much  larger  than  that  from  its  counterpart  of  glass-bead 
blasted  surface.  Note  that  the  force-deformation  curve  does  not  follow  the  initial  curve 
when  the  cylinder  is  re-loaded,  however,  the  slopes  of  the  initial  curve  and  the  reloaded 
curve  do  appear  to  be  the  same. 
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Figure  8.2:  A  sketch  of  experimental 


(Disk:  dia=  1  inch,  thickness  =  1/8  inch) 


Figure  8.3:  Top  view  of  a  surface  sandblasted  with  glass  beads  (Amplification:  x  100,  Scale 
Shown  as  the  white  bar  below  the  picture 
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Figure  8.4:  Profile  of  a  surface  sandblasted  with  glass  beads  -  a  side  view  (Amplification, 
and  scale  information  are  shown  in  the  dark  line  below  the  pictures) 
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Figure  8.5:  Variation  of  surface  asperities  traced  from  Figure  8.4. 
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Figure  8.6:  Top  view  of  a  surface  sandblasted  with  steel  beads  (Amplification:  x  100,  Scale. 
Shown  as  the  white  bar  below  the  picture). 
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Figure  8.7:  Profile  (at  two  locations)  of  a  surface  sandblasted  with  steel  beads  -  a  side  view 
(Amplification,  and  scale  information  are  shown  in  the  back  line  below  the  pictures) 


Figure  8.8:  Variation  of  surface  asperities  traced  from  Figure  8.7. 
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g  Figure  8.9:  A  steel  bead  wedged  in  the  asperities 
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8.1.4  Specimen  Preparation  and  Experimental  Arrangement  for  Determining 
the  Coefficient  of  Friction 

In  order  to  measure  the  values  of  the  coefficient  of  friction  for  the  surfaces  under  considera¬ 
tion,  a  simple  quasi-static  experiment  was  performed.  This  experiment  was  oriented  on  basic 
measurement  of  the  static  coefficient  of  friction,  as  identified  by  visible  inception  of  sliding. 
No  attempt  has  been  made  to  measure  initial  microslip,  change  of  normal  compliance  or 
other  more  intricate  phenomena  associated  with  frictional  sliding. 

The  disk  specimens  of  1/2  inch  thick  and  1  inch  diameter  were  cut  from  a  cold  rolled 
1020  steel  rod.  The  rod  is  from  the  same  stock  of  material  used  in  the  previous  test.  The 
material  has  an  elastic  modulus  E  =  30  x  106  psi,  a  yield  stress  approximately  42  x  103  psi. 
and  a  hardness  of  89  i?6- 

Both  surfaces  (top  and  bottom)  of  the  disk  specimen  were  prepared  in  a  same  manner 
as  the  previous  test.  The  disk  surfaces  were  first  lapped  with  a  milling  machine,  and  then 
sand-blasted  with  steel  beads  or  with  glass  beads  to  produce  an  artificial  surface  roughness. 
The  average  depth  of  surface  asperities  were  estimated  with  a  comparator  which  gave  a  value 
of  20  n  for  the  glass-bead-blasted  surface,  and  40  fx  for  the  steel- bead- blasted  surface  -  same 
as  those  reported  in  the  previous  tests. 

The  specimen  was  sandwiched  between  two  steel  anvil  blocks  which  have  a  hardened 
(approximately  32  Rc)  and  smooth  surface.  The  entire  assembly  was  then  installed  into 
the  test  apparatus.  A  normal  compressive  load  of  250  Lbf  was  first  applied  slowly  to  the 
assembly.  It  should  be  mentioned  that  due  to  the  smallness  of  the  surface  asperities,  we  were 
not  able  to  record  the  deformation  of  the  surface  asperities  during  loading.  After  completing 
the  normal  loading,  a  horizontal  load  is  slowly  applied  to  the  specimen  disk  until  the  disk 
begins  to  slip.  The  onset  of  slipping  was  observed  visually.  Two  tests  were  carried  out  for 
each  surface  roughness  and  the  results  are  summarized  in  the  next  section. 

Test  Results 


(a)  Glass-Beads-Blasted  Surface: 


Normal  load 
Horizontal  load 


=  250X6/ 

=  110X6/  (first  test) 

=  105X6/  (second  test) 
=  107.5X6/  (average) 


(8.37) 


Since  there  are  two  contact  surfaces  between  the  test  specimen  (roughened  on 
both  surface)  and  the  steel  anvil  blocks  (with  smooth  surfaces),  the  coefficient  of 
friction  is 
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(b)  Steel-Bead-Blasted  Surface 


n 


107.5 
2  x  250 


=  0.215 


Normal  load 
Horizontal  load 


The  coefficient  of  friction  is 


250L6/ 

106 Lbf  (first  test) 
lOOLbf  (second  test) 
103 Lbf  (average) 


(8.38) 


103 

2  x  250 


0.206 


It  should  be  mentioned  that  the  exact  magnitude  of  horizontal  load  at  the  onset  of  slipping 
is  difficult  to  ascertain  visually.  We  estimate  the  error  for  the  above  listed  horizontal  loads 
has  a  range  of  ±  5  Lbf. 


8.2  Numerical  Prediction  of  Interface  Contact 


We  used  the  new  microasperity  based  model  of  frictional  interfaces  to  predict  numerically 
the  force-deformation  relation  for  the  surfaces  studied  in  the  experiment.  The  computation 
was  performed  only  for  the  surface  sandblasted  with  glass  beads.  We  did  not  analyze  the 
surface  blasted  with  steel  beads,  because  it  was  very  rough  and  did  not  meet  at  least  two 
assumptions  of  our  model.  Those  assumptions  are:  small  deformations  of  asperities,  and 
non-interference  of  neighboring  asperities.  Moreover,  the  experimental  force- displacement 
curve  for  this  surface  indicates  that  the  error  of  measurements  on  unloading-loading  part 
was  rather  large  (see  Figure  8.11). 


8.2.1  Data  for  Numerical  Calculations 

To  calculate  statistical  characteristics  of  the  surface  we  used  scans  obtained  from  an  electronic 
microscope.  We  digitalized  the  side  view  of  the  sectioned  disk  specimen  in  aim  to  get  the 
surface  profile  as  a  function.  Graph  of  this  function  (height  versus  distance)  for  a  segment  of 
the  profile  is  shown  in  Figure  8.12.  The  figure  also  shows  the  actual  scan  of  the  corresponding 
segment  of  the  real  profile. 

Using  the  sampling  method,  we  obtained  the  following  statistical  characteristics  of  the 
surface: 
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Figure  8.12:  Digitized  and  original  profile  of  the  specimen  used  in  the  experiment. 


140 


Mean  height: 

112 

fim, 

Mean  Slope: 

0.000729 

o 

Ci 

II 

Mean  Curvature: 

-0.000288 

Atm'1, 

Standard  deviation  of  height: 

4.91 

Hm, 

Standard  deviation  of  slope: 

0.160 

Standard  deviation  of  curvature: 

0.0227 

Surface  peak  density: 

0.000508 

peaks/ fim2 

Wavelength  spectrum  parameter: 

0.304. 

Another  microscopic  scan  of  the  specimen  (Figure  8.3)  indicates  that  the  surface  is 
isotropic  so  that  the  above  coefficients  fully  characterize  the  shape  of  the  surface. 

To  be  able  to  perform  mechanical  analysis  of  asperities  we  also  need  to  know  material 
constants  of  the  specimen  material  (cold  rolled  SAE  1020  steel).  Bodner-Partom  constitutive 
constants  for  steel  are  given  in  reference  [10].  These  generic  constants  were  tuned  for  our 
particular  sample,  using  the  results  of  the  tension  test  (Figure  8.13).  The  final  values  of  the 
constants  used  for  numerical  calculations  are  listed  below. 


E  =  209  GPa 
mi  =  0.05  MPa 
n  =  2.28 
v  =  0.3 
Z0  =  600 MPa 
D0  =  10000s'1 
Zi  =  1050  MPa 

These  constants  provide  good  correlation  of  numerical  and  experimental  results  for  the 
tension  test  -  see  Figure  8.13. 

8.2.2  Modeling  of  Surface  Loading 

As  a  first  step  of  the  homogenization  procedure,  the  above  surface  statistics  was  used  to 
determine  effective  support  of  joint  probability  density  of  asperity  peaks,  in  particular  the 
range  of  surface  peak  heights  and  curvatures.  This  range  is  defined  by: 
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Figure  8.13:  Comparison  of  tensile  stress-strain  relationships  obtained  experimentally  and 
numerically. 
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Figure  8.14:  Adaptively  refined  mesh  for  an  axisymmetric  asperity  (deformed  shape  is 
shown). 


—12 nm  <  zp  <  16 nm  (8.39) 

0.004/im-1  <  kp  <  0.08/xm~1  (8.40) 

Within  this  range,  several  sample  asperities  under  normal  load  were  modeled  by  the  finite 
element  method.  Because  the  experimental  surface  was  very  rough  and  no  effort  has  been 
made  to  make  it  extremely  clean,  adhesion  effects  were  neglected.  It  was  assumed  that  the 
deformation  of  the  asperity  was  elasto- viscoplastic. 

Figure  8.14  shows  a  deformed  mesh  which  was  used  to  discretize  half  of  the  asperity 
cross  section  (the  problem  is  axisymmetric).  The  rate  of  load  application,  controlled  by  the 
normal  velocity,  was  equal  to  0.25  /im/s,  and  corresponded  roughly  to  the  one  used  in  the 
experimental  setup. 

About  10,000  time  steps  were  necessary  to  complete  the  analysis.  At  each  of  those  points 
the  values  of  approach  contact  force  and  area  of  contact  were  printed  out.  From  those  10,000 
discrete  values  about  40  were  chosen  as  data  for  statistical  homogenization.  This  was  because 
not  all  of  those  10,000  values  were  computed  with  the  same  accuracy.  It  can  be  explained 


area  [/mi2] 


Figure  8.15:  Area  of  contact  as  a  function  of  approach  obtained  directly  from  the  numerical 
analysis  (about  10,000  data  points). 

on  the  example  of  area  of  contact.  As  a  result  of  numerical  analysis  we  obtain  the  area 
in  terms  of  approach  as  a  piecewise  constant  function  (see  Figure  8.15)  while  it  is  obvious 
that  this  relationship  is  at  least  continuous  for  a  smooth  asperity.  The  discontinuity  of  the 
numerical  result  is  caused  by  the  fact  that  the  contact  condition  is  examined  at  integration 
points.  Therefore,  the  contact  area  increases  incrementally,  as  new  integration  points  join 
the  contact  zone.  For  further  analysis  the  centroids  of  the  horizontal  segments  on  the  graph 
of  area  versus  approach  were  chosen  as  basic  data  points.  Final  graphs  of  contact  force  and 
area  in  terms  of  approach  axe  shown  in  Figures  8.15  and  8.17,  respectively. 
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Figure  8.17:  Final  graph  of  contact  force  as  a  function  of  asperity  deflection. 

The  results  of  mechanical  analysis  shown  in  Figures  8.15  through  8.17  were  obtained  for 
one  representative  asperity  with  width  to  height  ratio  w/h  =  13.4.  For  exactly  cosinusoidal 
profile  this  w/h  ratio  gives  standard  deviation  equal  to  standard  deviation  of  slope  of  the 
real  surface.  The  wavelength  spectrum  parameter  of  the  analyzed  surface  indicates  that  also 
asperities  with  other  w/h  could  have  an  influence  on  the  expected  values  of  global  force  and 
area  of  contact.  However,  for  the  sake  of  simplicity  and  efficiency  of  the  computation  we 
limited  our  numerical  analysis  to  the  asperities  with  the  most  representative  width  to  height 
ratio.  Note  that,  due  to  absence  of  rate-dependent  effects,  the  results  obtained  for  the  basic 
asperity  (peak  curvature  equal  to  0.042  /im'1)  are  also  useful  for  asperities  with  different 
peak  curvatures.  The  contact  force  and  area  for  those  asperities  can  be  computed  by  simple 
rescaling  of  the  existing  values,  because  both  force  and  area  of  contact  axe  proportional  to 
squared  dimensions  of  an  asperity. 

The  results  of  finite  element  asperity  analysis,  in  particular  contact  force  and  read  contact 
area,  were  used  by  statistical  homogenization  package,  to  produce  force-approach  and  area- 
pressure  curves  shown  in  Figures  8.18  and  8.19.  Note  that,  at  this  stage,  we  focus  only  on  the 
loading  part  of  the  curve.  The  numerical  predictions  are  compared  with  the  experimental 
results.  Note  that  the  original  experimental  curve  was  rescaled  to  represent  only  one  surface 
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Figure  8.18:  Numerical  and  experimental  approach-pressure  curves  for  the  interface. 

pair  and  normal  pressure  rather  than  total  force.  In  the  spirit  on  remark  from  Section  6.1, 
the  zero  point  on  the  theoretical  “approach”  axis  was  chosen  to  match  experimental  load 
and  surface  stiffness  at  the  initial  contact  of  the  surfaces. 
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A  comparison  of  experimental  and  numerical  results  shows  a  very  good  agreement  up  to 
a  load  of  about  800  N/^m2  (1000  psi).  For  higher  loading  an  increasing  discrepancy  of  the 
experimental  and  numerical  curves  can  be  observed  (see  Figure  8.20).  There  can  be  several 
reasons  of  that  and  some  of  them  are  discussed  below. 

1.  To  predict  more  precisely  behavior  of  the  specimen,  the  elastic  deformation  of  the 
bulk  material  has  to  be  considered.  A  graph  of  experimental  normal  stress  reduced 
by  bulk  deformation  is  shown  in  Figure  8.21.  Apparently,  for  higher  loading  bulk 
deformation  is  significant  and  has  to  be  taken  into  account. 

2.  Another  reason  of  discrepancy  between  numerical  and  experimental  results  can  be 
additional  deformation  of  the  apparatus  and  stacked  samples.  Assuming  pro¬ 
portionality  of  that  deformation  to  the  loading  one  can  tentatively  modify  experimental 
data.  Figure  8.22  shows  original  and  modified  experimental  results.  It  was  assumed 
that  the  experimental  measurement  has  systematic  error  which  depends  linearly  on 
loading  and  causes  that  the  sensor  indicates  displacement  exaggerated  by  1.5  /im  for 
the  highest  loading.  In  Figure  8.23,  the  modified  experimental  data  are  compared  with 
numerical  results.  A  good  agreement  of  both  curves  can  be  observed  for  wide  range  of 
loading,  which  confirms  strong  possiblity  of  systematic  settling  in  the  apparatus  and 
in  the  samples. 

3.  Error  of  numerical  results  increases  with  increasing  loading  because  at  higher  loads 
more  asperities  have  large  deformations,  so  that  more  asperities  are  beyond  the  assump¬ 
tions  of  our  theory  (see  Figure  8.24). 

4.  Another  reason  of  the  discrepancy  for  loads  higher  than  800  N/cm2  can  be  loss  of 
accuracy  during  experiment  caused  by  unloading.  Complete  unloading  might  have 
been  accompanied  by  a  microdisplacement  in  a  horizontal  direction  and  then  reloading 
would  deform  different  asperities  than  the  original  loading. 
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Figure  8.20:  Numerical  and  experimental  contact  force  in  terms  of  approach  for  wider  range 
of  loadings. 


n.stress  [107  Pa] 


Figure  8.21:  Experimental  stress  versus  approach  curves:  (a)  original  and  (b)  corrected  for 
bulk  deformation  of  a  sample. 
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Figure  8.22:  Experimental  pressure-approach  curves:  (a)  original  and  (b)  corrected  for  ad¬ 
ditional  settling  of  the  setup. 
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8.3  Studies  of  Unloading 


An  additional  study  that  was  performed  for  the  experimental  surface  was  focused  on  modeling 
of  the  static  unloading  path.  By  static  unloading  we  understand  the  situation,  wherein 
the  normal  load  is  reduced  to  zero  or  almost  zero  without  relative  sliding  on  the  interface. 
This  means  that  the  asperity  tips,  which  became  “conformal”  during  loading  due  to  plastic 
deformation,  remain  conformal  and  aligned.  This  situation  is  illustrated  in  figure  8.25. 

Importantly,  modeling  of  the  unloading  path  required  additional  extension  of  the  homog¬ 
enization  package  to  monitor  the  displacement  and  load  corresponding  to  the  beginning  of 
surface  unloading.  Moreover,  even  for  single  unloading  of  the  interface,  the  results  at  the 
asperity  level  must  consider  a  whole  family  of  unloading  curves,  starting  at  different  normal 
deflection.  This  is  a  consequence  of  random  asperity  height;  at  the  beginning  of  surface 
unloading  each  individual  asperity  is  “caught”  at  a  different  stage  of  deformation. 

For  practical  modeling,  it  would  not  be  feasible  to  consider  an  infinite  number  of  unload¬ 
ing  paths  at  the  asperity  level.  Instead,  we  adopted  the  following  procedure: 

(a)  analyze  in  detail  a  few  (three  to  four)  unloading  paths  for  the  asperity,  starting 
at  different  deflection  levels,  and 

(b)  use  approximation  techniques  to  represent  unloading  paths  which  start  at  the 
intermediate  load  level. 

This  approach  is  illustrated  in  figure  8.26.  Importantly,  the  finite  element  analysis  of  the 
unloading  path  was  performed  using  exactly  the  same  methods  as  for  the  loading  curve  - 
the  unified  viscoplastic  constitutive  theories  applied  here  do  not  require  special  treatment  of 
elasto-plastic  unloading.  The  approximation  method  used  for  the  construction  of  intermedi¬ 
ate  unloading  paths  was  based  on  the  blending  function  formula  [  ]. 

The  formula  is  a  mapping  of  a  unit  square  into  a  domain  which  has  four  curvilinear  sides. 
Coordinates  of  the  square  (s,  t  e  [0,1])  can  be  treated  as  parameters.  They  parametrize  sides 
of  the  curvilinear  domain.  Let  the  parametric  equations  of  these  curves  71  =  AB,  72  =  BC, 
73  =  DC,  74  =  AD  (Figure  8.26)  be  given  generally  in  the  form 

7 k  =  {(*,  y)&2  :  x  =  xk(p),  y  =  yk(p),  pe{ 0, 1]}  ,  k  -  1, 2, 3, 4 

where  in  our  case: 


x  are  approaches, 

y  are  contact  forces, 
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Figure  8.25:  Static  loading,  unloading  and  reloading  of  the  interface 


s  for  k  =  1,3, 
t  for  k  =  2, 4. 


Using  the  above  notations  the  blended  functions  provide  mapping  of  parameters  s,  t  into 
x  and  y  coordinates  in  the  following  way 


z(s,  t)  =  xi(s)(l  -  t)  4-  x2(t)  ■  s  +  x3(s)  •  t  + 

+x4(t)  •  (1  -  s)  -  xA{\  -  t)(  1  -  s)  -  xB  ■  (1  -  t)  •  s 
—Xc  ■  t  ■  3  —  Xd  •  t  ■  (1  —  s) 

y(s,t)  =  t/i(s)(l  -  t)  +  y2(t)  •  s  +  y3(s)  ■  t  + 

+y4{t)  •  (1  -  s)  -  y^(l  -  t)(  1  -  s)  -  yB  ■  {l  -  t)  ■  s 
-yc  •  t  ■  s  -  yD  ■  t  •  (1  -  s) 

where  (xA,yA),(xB,yB)Axc,yc)AxD,yD)  are  coordinates  of  corners  A,  B,  C,  D  correspond¬ 
ingly.  Figure  8.27  shows  a  family  of  unloading  curves  generated  by  the  blending  functions. 
Note  that  the  blended  curves  can  be  generated  also  outside  domain  ABCD.  Curve  70  was  not 
generated  but  computed  by  solution  of  the  unloading  problem  for  the  asperity.  Comparison 
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Figure  8.26:  Contact  force  for  an  asperity  in  terms  of  deformation  during  loading  and  un¬ 
loading  at  three  different  levels. 
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of  this  curve  with  blended  functions  indicates  reasonably  good  correlation.  This  correlation 
improves  with  increasing  number  of  calculated  unloading  curves. 

After  application  of  the  statistical  homogenization,  the  loading-unloading  curves  were 
obtained  for  the  surface,  and  are  shown  in  figure  8.28.  Clearly,  the  theoretical  unloading 
path  is  steeper  than  the  one  observed  in  the  experiments.  At  present  time  we  do  not  have 
a  precise  explanation  of  this  fact.  However,  we  believe,  that  the  asperity- based  theoretical 
curve  better  represents  the  actual  behavior  of  the  interface.  This  is  supported  by  the  following 
arguments: 

1.  For  the  loading  section  of  the  curve,  random  ^perity  height  has  a  significant  smoothing 
effect,  because  different  asperities  are  contacted  at  different  loading  stages.  However, 
the  unloading  path  is  only  supported  by  elastic  “rebound”  of  all  the  deformed  asperities 
-  see  figure  7.26.  Therefore,  the  unloading  deflection  cannot  essentially  be  higher  than 
the  maximum  elastic  rebound  of  a  single  asperity.  Since  the  unloading  path  for  a 
single  asperity  is  rather  steep  (see  figure  8.26),  the  elastic  rebound,  even  for  the  largest 
asperities,  does  not  exceed  1.0  /j.. 

2.  The  additional  normal  compliance,  present  in  the  experimental  setup,  can  have  a  much 
more  significant  distorting  effect  on  the  (very  steep)  unloading  path  than  on  the  load¬ 
ing  path.  Indeed,  presence  of  additional  deflection  of  about  1.5  fi  would  modify  the 
theoretical  unloading  curve  to  match  the  experimental  one.  It  is  conceivable  that  such 
deflection  can  be  present  due  to  departure  of  the  surfaces  from  ideal  flat,  and  the 
resulting  “plate  bending”  effect. 

It  is  also  possible,  that  there  may  be  some  infinitesimal  tangential  sliding  present  on  the 
surface  during  unloading  -  this  would  spoil  the  alignment  of  deformed  asperities  shown  in 
figure  8.25  and  cause  additional  separation  of  the  surfaces. 

The  above  arguments  are  supported  by  the  experimental  work  of  Connoly,  Shoffield  and 
Thornley  [30],  wherein  experimental  unloading  curves  were  much  steeper  than  these  obtained 
in  the  experiment  reported  here. 

In  view  of  these  remarks,  the  question  of  interface  unloading  will  require  some  additional, 
extremely  precise  experimental  studies.  Furthermore,  the  phenomenon  of  unloading  during 
tangential  sliding  will  need  a  dedicated  study.  The  unloading  curve  in  this  case  will,  due  to 
mechanical  interaction  of  the  sliding  asperities,  differ  significantly  from  the  static  case. 


158 


6.00 


8.00  10.00 


approach [n 


curves  obtained  numerically  and  experimentally. 


9  Studies  of  Friction 


In  this  section,  we  present  some  introductory  studies  of  frictional  characteristics  of  the  in¬ 
terface.  In  particular,  we  aim  to: 

•  estimate  the  value  of  the  static  coefficient  of  friction,  and 

•  study  frictioncal  deformation  of  microasperities. 

9.1  Static  Coefficient  of  Friction 

For  the  sample  surfaces  studied  in  previous  sections,  introductory  studies  of  the  coefficient 
of  friction  were  performed.  These  studies  were  based  on  the  following  two  assumptions: 

1.  sliding  resistance  between  similar  metallic  surfaces  is  caused  primarily  by  formation 
and  shearing  of  junctions  between  the  contacting  asperities, 

2.  the  value  of  tangential  force  necessary  to  shear  a  junction  is  defined  as: 

T  =  t,Ac  (9-41) 

where  r,  is  shear  resistance  of  the  material  and  Ac  is  the  real  contact  area  between  the 
asperities. 

In  the  first  approximation  the  shear  strength  of  a  junction  can  be  considered  to  be  indepen¬ 
dent  of  normal  pressure  and  equal  to  the  shear  strength  of  the  bulk  material. 

The  formula  (9.41)  can  be  easily  used  to  calculate  the  coefficient  of  friction  on  the  interface 
level: 


fi  =  TaAc/F 

where  Ac  is  a  real  contact  area  per  unit  nominal  area  and  F  is  the  normal  force  per  unit 
nominal  area. 

Smooth  Surface 

The  above  formula,  combined  with  pressure-area  curves  presented  for  the  smooth  alu¬ 
minum  surface  in  Figure  7.9,  produces  the  values  of  the  coefficient  of  friction  shown  for 
increasing  normal  load  in  Figure  9.1.  For  the  case  with  adhesion  (clean  smooth  surfaces  in 
vacuum),  the  value  of  the  coefficient  of  friction  is  very  high  at  light  loads  -  this  is  caused  by 
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Figure  9.1:  Smooth  aluminum  surface:  coefficient  of  friction  for  increasing  loads 

additional  normal  attraction  due  to  adhesion.  In  the  absence  of  adhesion  the  coefficient  of 
friction  is  almost  constant  over  a  wide  range  of  loads,  with  a  value  slightly  below  0.5.  Note 
that  this  value  is  essentially  overestimated  -  junction  strength  of  oxygenated  aluminum  sur¬ 
faces  is  lower  than  the  shear  resistance  of  the  bulk  material,  assumed  here.  A  more  detailed 
studies  of  this  phenomenon,  including  frictional  microdisplacements,  are  currently  underway. 

Rough  Surface 

Some  interesting  observations  can  be  made  when  the  above  analysis  is  applied  to  the 
rough  surface,  studied  in  Section  7.4.  The  values  of  the  coefficient  of  friction  calculated 
for  this  surface  using  elastic  and  elasto-plastic  models  are  compared  in  Figure  9.2.  Clearly, 
a  purely  elastic  approach  underestimates  the  value  of  the  coefficient  of  friction.  This  is  a 
result  of  an  underestimated  contact  area  predicted  by  the  elastic  solution  on  the  asperity 
level.  The  elasto-plastic  approach  produces  more  realistic  values  of  the  coefficient  of  friction, 
of  the  order  of  0.3. 

Experimental  Surface 

The  above  simplified  model  was  used  to  estimate  the  coefficient  of  friction  for  the  rough 
steel  surface  studied  experimentally  and  numerically  in  Section  8.  Figure  9.3  shows  the 
calculated  value  of  the  coefficient  of  friction  for  increasing  normal  load  (approach).  In  the 
absence  of  relevant  data,  it  was  assumed,  that  the  shear  resistance  of  asperity  junctions  is 
equal  to  shear  strength  of  the  bulk  material.  For  the  actual  normal  pressure  tested  in  the 
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experiment,  the  estimated  value  of  the  coefficient  of  friction  due  to  junction  shearing  was 
equal  to  0.093,  while  the  measured  value  was  0.215. 

However,  for  this  rather  rough  surface  there  exists  an  additional  contribution  to  friction 
due  to  interlocking  of  asperities.  While  for  relatively  smooth  surfaces  this  component  can 
usually  be  neglected,  for  this  rough  surface  this  is  not  the  case.  A  very  approximate  estimate 
of  the  interlocking  component  of  friction,  based  on  relative  motion  of  inclined  surfaces  [32,89], 
leads  to  a  conclusion  that  the  additional  contribution  to  the  coefficient  of  friction  is  equal  to 
the  average  slope  of  the  asperities.  For  the  surface  under  consideration  this  is  equal  to  0.16, 
which  leads  to  the  final  value  of  the  coefficient  of  friction  equal  to  0.253,  as  compared  with 
experimental  value  of  0.215. 

This  is  a  quite  good  correlation  of  numerical  and  experimental  results,  especially  that, 
due  to  a  lack  of  relevant  data,  the  shear  resistance  of  asperity  junctions  was  assumed  to  be 
independent  of  the  normal  pressure  (which  is  not  exactly  the  case  in  practice). 

9.2  Studies  of  Frictional  Sliding 

The  simplified  approach  discussed  in  the  previous  subsection  can  only  be  used  for  general 
estimates  of  sliding  resistance  of  the  interface.  It  cannot  represent  such  intricate  phenomena 
as: 


•  tangential  micro-displacements  before  inception  of  macro-scale  sliding, 

•  damage  of  asperities  by  shearing  below  the  contact  junction,  rather  than  at  the  junction, 

•  normal  “settling”  of  the  junction  due  to  combined  normal  and  tangential  stress, 

•  ploughing  component  of  friction,  etc. 

Studies  and  understanding  of  the  above  phenomena  can  only  be  accomplished  by  mod¬ 
eling  of  the  actual  sliding  process,  wherein  two  asperities  in  contact  are  subject  to  tangen¬ 
tial  mot’on.  Numerical  modeling  of  these  phenomena  is  computationally  rather  expensive, 
because  mil  three-dimensional  asperity  models  need  to  be  used,  combined  with  extensive 
nonelastic  deformation,  complex  stress  states  and  sliding  resistance  of  the  metallic  junction. 

In  this  section,  some  introductory  studies  of  tangential  sliding  are  presented. 

We  analyzed  a  representative  cosinusoidal  asperity  with  height  equal  to  2  microns  and 
diameter  at  the  base  equal  to  10  microns.  The  asperity  was  made  of  the  steel  studied  in 
the  experiment,  with  full  elasto-viscoplastic  representation.  The  contacting  flat  was  first 
subjected  to  normal  motion  (0.125  microns),  followed  by  tangential  sliding  of  1.0  micron. 
The  first  step  was  performed  in  time  0-1800  s.  The  tangential  loading  step  was  applied  in  time 
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Figure  9.4:  Evolution  of  normal  and  tangential  forces  for  the  asperity 


1800-4000  s.  In  order  to  model  shear  resistance  of  the  metallic  junction,  large  coefficient  of 
friction  (20)  was  assumed.  The  main  objective  here  was  to  identify,  whether  the  maximum 
deformation  would  occurr  on  the  surfaace  (shearing  of  the  junction)  or  below  the  surface 
(shearing  ot  the  asperity  tip). 

The  asperity  was  discretizated  by  FEM  mesh  which,  after  h-adaptation,  had  666  degrees 
of  freedom  Introductory  results  of  this  study  are  shown  in  Figs.  9.4  and  9.5.  In  particular, 
figure  9.4  shows  evolution  of  normal  and  tangential  forces  during  the  loading  process.  Figure 
9.5  shows  deformation  of  the  asperity  at  the  end  of  the  analysis  (at  time  4000s). 

During  the  analysis,  no  signs  of  destruction  of  the  asperity  were  observed.  In  fact,  the 
average  tangential  traction  on  the  contact  surface  due  to  tangential  motion  of  the  flat  was 
approximately  64MPa,  which  is  well  below  the  shear  strength  of  steel.  This  means,  that 
representing  shear  strength  of  the  junction  by  high  coefficient  of  friction  was  not  sufficient. 
Instead,  exact  modeling  of  shear  resistance  on  the  junction  surface  is  necessary.  This  required 
additional  modification  of  the  finite  element  code  and  study  of  asperity  sliding.  These  studies 
are  currently  underway  and  will  be  presented  in  the  forthcoming  publications.. 
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Figure  9.5:  Deformed  configuration  and  contours  of  horizontal  displacement  in  the  asperity 
sliding  study. 
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10  Towards  Application  of  Asperity-Based  Models  in 
Modeling  of  Dynamic  Friction 

The  micro-asperity-based  models  developed  in  this  work  produce  the  values  of  normal  pres¬ 
sure,  real  contact  area,  sliding  resistance,  etc.,  as  a  function  of  approach  of  the  two  surfaces. 
These  functions  are  given  in  the  form  of  series  of  discrete  points  or  graphs. 

However,  application  of  these  constitutive  models  in  the  analysis  of  dynamic  friction  re¬ 
quires  an  analytic  (and  differentiable)  formulas.  Towards  this  end,  some  effort  was  dedicated 
to  development  of  automated  tools  that  would  produce  analytical  formulas  representing  the 
interface  constitutive  laws. 

A  basic  approach  used  is  as  follows: 

1.  micro-asperity  based  model  is  given  as  a  series  of  discrete  points,  for  example  a  sequence 
of  approach-pressure  values  (Ak,Pk)  k  ~  1,M 

2.  the  user  chooses  the  formula  to  be  used,  say  the  Oden-Martins  law  or  logarithmic 
compliance  curve, 

3.  the  coefficients  of  the  constitutive  formula  are  determined  automatically  through  the 
minimization  procedure  described  below. 

The  above  procedure  will  be  illustrated  on  the  example  of  the  fitting  of  the  asperity- 
based  approach-pressure  relationship  into  the  Oden-Martins  interface  model  [71].  The  Oden- 
Martins  normal  compliance  law  can  be  written  in  the  following  generalized  form: 

p  =  c(^r  (10.42) 

ao 

where 

p  —  nominal  pressure, 

a  —  approach, 

c,g,m  —  unknown  material  constants, 

ao  =  1  pm  (dimensionality  factor). 

Note  that  ao  was  introduced  to  assure  proper  dimensionality,  and  g  is  a  translation  (initial 
gap)  which  appropriately  locates  the  zero  point  for  the  approach  axis.  The  constants  in  the 
above  formula  were  determined  in  such  a  way  which  provides  the  best  fitting  of  the  function 


167 


pressure  [104  Pa] 


Figure  10.1:  Fitting  of  the  Oden  -  Martins  formula  to  the  sequence  of  asperity- based  data 
points  for  the  steel  surface  studied  in  Section  8 

10.42  to  the  sequence  of  M  given  points  (a.k,Pk)-  The  best  fitting  means  here  that  a  function 
B,  which  is  square  of  a  norm  of  differnce  between  discrete  and  continuous  laws,  attains  its 
infimum. 

Function  B  has  the  following  form: 


M 

B(c,g,m )  =  Y,lP(c'9,™,ak)  -  Pk]2 


fc=  i 


(10.43) 


where 


p(c,g,m,a)  —  nominal  pressure  given  by  formula  10.42 
( ak,pk )  —  k-th  discrete  point. 

For  the  load-approach  curve  representing  our  experimental  surface  (see  Section  8),  we 
obtained  g  =  0 pm,  c  =  1 5.2N/cm2,  and  m  =  2.25.  Figure  10.1  shows  pressure  -  approach 
relationships  for  both  discrete  and  the  best  fitting  continuous  forms. 

A  similiar  procedure  was  applied  to  the  hypothetical  very  smooth  surface  described  in 
section  7.  We  obtained  g  =  0.0025pm,  c  =  3.17  x  10 9N/cm2,  and  m  =  5.17.  The  graphs  of 
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pressure  [107  Pa] 


Figure  10.2:  Fitting  of  the  Oden  -  Martins  formula  to  the  sequence  of  discrete  points  for  the 
very  smooth  surface 

aseperity-based  and  analytical  pressure-approach  relationships  are  shown  in  figure  10.2. 

These  examples  indicate,  that  Oden- Martins  law  provides  a  good  representation  of  typical 
constitutive  models  of  metallic  interfaces.  Note  that  the  approach  presented  here  is  applicable 
to  other  formulations,  such  as  exponential  normal  compliance  law. 
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11  Conclusions 


This  report  presented  a  detailed  theoretical  background,  numerical  algorithm,  and  practical 
results  of  the  development  of  new  asperity- based  constitutive  models  of  frictional  interfaces. 
While  the  main  idea  stems  from  previous  works  on  asperity-based  approaches,  the  application 
of  the  adaptive  finite  element  method  to  the  modeling  of  nonlinear  asperity  deformation 
brought  these  models  to  a  level  of  practical  application  to  the  modeling  of  real  engineering 
surfaces.  This  is  confirmed  by  successful  comparisons  of  the  asperity- based  predictions  with 
the  results  of  carefully  designed  verification  experiments. 

The  new  models  developed  in  this  work  will  find  application  in  many  diverse  aspects  of 
the  modeling  of  friction,  such  as: 

•  simulation  and  control  of  friction- induced  squeaks,  stick-slip  motion,  chatter,  and  other 
unstable  phenomena, 

•  precise  modeling  of  tribological  surfaces,  such  as  rollers,  bearings,  etc. 

•  modeling  the  conductivity  of  thermal  and  electrical  connections,  including  microelec¬ 
tronic  devices  and  semiconductors, 

•  understanding  and  modeling  of  surface  damage  and  wear  mechanisms, 

•  and  many  others. 

Indeed,  the  results  of  this  research  project  are  already  finding  their  way  into  practical 
applications  in  several  of  the  above  areas. 

It  is  important  to  note,  however,  that  there  remain  many  open  questions  and  challenges 
in  this  topic,  including: 

•  extensions  to  hyperelastic  and  brittle  materials, 

•  extensive  studies  of  surface  sliding  and  various  components  of  frictional  resistance 
(shearing,  interlocking,  and  ploughing), 

•  studies  of  mechanisms  of  surface  damage  and  wear, 

•  extensions  to  hydrodynamics  lubrication, 

•  dynamic  loading  of  the  surface  and  high-velocity  sliding,  etc. 

These  issues  need  to  be  addressed  in  the  future  to  provide  a  full  understanding  and 
control  of  the  complex  phenomena  occurring  in  frictional  interfaces. 
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